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1.0 INTRODUCTION 


This report contains the documentation and user's guide for the Analytical Satel- 
lite Orbit Predictor (ASOP) computer program. ASOP is based on mathematical 
methods that represent a new state-of-the-art for rapid orbit computation tech- 
niques. The theoretical development of these methods has been carried out dur- 
ing the past few years, and they are now in the form of an operational computer 
program. ASOP is intended to be used for computation of near -Earth orbits 
including those of. the Shuttle/Orbiter and its payloads. 

Orbit computation methods can usually be classified as: 

a. Numerical methods - The calculations are carried out in a step-by-step man- 
ner. High precision is possible, but computer runtime can be exessive. 

b. Analytical methods - The- calculations are carried out in one step regardless 
of the prediction interval. Therefore, these methods have extremely fast 
computation times . 

In the past, analytical methods have not been widely used because they were less 
accurate and required much more computer coding than numerical methods. The 
Poincare -Similar elements (PS<|>) used in ASOP overcomes these disadvantages. It 
is possible to compute near-Earth orbits to within an accuracy of a few meters. 
Recursive equations are used instead of complicated formulas. Execution time of 
ASOP is on the order of a few milliseconds. 

The theoretical foundation for the mathematical techniques used in ASOP were 
developed by Dr. Gerhard R. Scheifele (refs. 1 and 2). Scheifele developed the 
Delaunay-Similar elements (DS<j>) based on the true anomaly and used the elements 
to solve the J 2 perturbation problem (ref. 3). Later, Scheifele developed 
the Poincare-Similar elements (PS<j>) , which contain no singularities for zero 
eccentricity and zero inclination. 

In reference 4, Mueller describes the relationship between the PS<j> elements and 
the Cartesian coordinates and establishes the PS4» perturbed equations of motion. 
These elements and the associated equations of motion were used by Bond (ref. 5) 
to develop a nonsingular analytical solution to the J 2 perturbation problem. 

In 1977 the analytical solution was expanded by Scheifele (ref. 6) and Mueller 
(ref. 7) to include the perturbations due to atmospheric drag. Later, Mueller 
(ref. 8) developed the higher order zonal geopotential terms that have been 
implemented and documented by Wang (ref. 9). The development of time-dependent 
(tesseral) geopotential perturbation theory has been completed and is described 
by Mueller in reference 10. All of the current analytical theory has been 
implemented into ASOP including the zonal long period geopotential terms. 

The ASOP program has been designed in two versions: (1) a stand-alone version 

that can be used interactively to obtain immediate results to a specific problem 
and (2) a user-subroutine package that can be incorporated into other software 
systems. Both versions were designed to be small and to execute quickly. 
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Both versions of the ASOP program have been written in UNIVAC standard FORTRAN- 
V and are available to the public on file NUMEG under the qualifier FM 6 -NO 8569 
on th'e UNIVAC 1110-Exec 8 system. This document is intended to instruct the 
user in the operation of the ASOP program on this machine and to document the in- 
dividual ASOP subroutines. 
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2.0 USER'S GUIDE 

This section is intended to give the user all the information necessary to 
operate the ASOP programs. Because the program is designed to operate in two 
modes (stand-alone and subroutine package), each mode of operation is described 
separately. 

The first part of this section (sec. 2.1) will describe the general input pa- 
rameters and options available when using the stand-alone ASOP program. Also 
described in this section are the standard default values and the typical com- 
mands needed to execute the program in the demand or the batch mode . Finally 
some sample output is given to help the user if modifications are to be made to 
the system. 

Section 2.2 will deal with the ASOP subroutine package. This section will de- 
scribe the necessary modules that are used within the package as well as any 
interface requirement that, the user must be aware of if he is to include this 
package in his own software . The input to and output from the ASOP subroutine 
are also fully described in this section, as are the subroutine’s default values. 


2.1 INDEPENDENT PROGRAM 

The ASOP program was designed as an interactive program capable of giving the 
user fast , accurate answers to Shuttle- type orbit problems . The program may 
also be run in a batch environment if a large number of cases are to be 
investigated. 

There are two basic methods used to control the operation of the ASOP program: 
flags and direct-user interaction. The flags are used to indicate the type of 
data being entered and to select certain options within the program. Direct in- 
teraction allows the user to check the input data to ensure their accuracy 
before continuing. 

Primary data input to the ASOP program is accomplished using the NAMELIST 
'$INPUT' . The necessary input variables and user responses to program questions 
are described in section 2.1.1, and the program default values are described in 
section 2.1.2. Section 2.1.3 explains the printed output generated by the ASOP 
program. Finally, sections 2.1.4 and 2.1.5 describe the commands and user re- 
sponse required to run the ASOP program and give an example of the resulting 
output . 


2.1.1 Input Description 

The NAMELIST is the primary method of getting data into the ASOP program. How- 
ever, during normal operation, the user is expected to interact with the program 
by supplying additional information. After starting the ASOP program (sec. 
2.1.4), it will ask for the NAMELIST data with the statement 


INPUT DATA USING NAMELIST '$INPUT‘ 
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At this point, the user has three options: 

a. To enter the NAMELIST data directly from the keyboard. 

b. To add a data file or element containing the NAMELIST information using the 
§ADD command (ref. 11). 

c. To enter @E0F to terminate program execution. 

If option C is selected, the program will respond with 

**NORMAL PROGRAM TERMINATION** 


and the program will stop.. If option A or B is selected the program will print 
out all the NAMELIST variables and their associated values (including default), 
as well as the initial conditions of the problem. The message 


ENTER: X = EXECUTE; S = STOP; C r CHANGE INPUT 


should then appear. Here, the user should check the input data and enter the 
necessary letter. If an X is entered, the program will continue the execution 
as directed by the input. When the input stop condition is satisfied, the pro- 
gram will again ask for data input as described earlier. The series of instruc- 
tions can be repeated as often as necessary. 

If an S is entered in response to the message, then the program will respond 
with **N0RMAL PROGRAM TERMINATION** and program execution will terminate. 

Should a C, or any other alphanumeric character, be entered at this point, 
then the message 


**CHANGE DATA USING THE NAMELIST '$INPUT'** 


will be displayed. The user can then reenter those values that are incorrect. 
Once the necessary corrections have been made, the program will again display 
the input data, the initial conditions, and the message 


ENTER: X = EXECUTE: S = STOP: C = CHANGE INPUT 


Table I describes the input variables that may be used in the NA MEL IST INPUT' . 
Whether keying in the information or creating a data element, a H$INPUT must be 
entered first where the K represents one or more spaces. • Each variable 
entered must be preceded by one or more spaces, and if more than one variable is 
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to appear on a line, they 
nate the NAMELIST input a 
See reference 12, pages 6 
statement. 


must be separated by a blank or a comma ( , ) . To termi- 
K$END or must be the last item entered, 

through 13, for a complete description of a NAMELIST 
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TABLE I.- NAMELIST INPUT VARIABLES 


Variable Type 

EL DP 


IEL 


STOP 


ISTOP 


DP 


PRINT 


DP 


No. of inputs 
required 


Description and 
available option s 


Must be supplied by the user; may be 
given in the following forms as deter- 


mined by 

the flag 

IEL 


EL (1) 

= 

X 

or 

a 

or 

»a b 

EL (2) 

= 

Y 

or 

e 

or 

hp b 

EL (3) 

= 

Z 

or 

i 



EL (i|) 

= 

• 

X 

or 

03 



EL (5) 

= 

# 

Y 

or 

a 



EL (6) 


Z 

or 

M 




All angular input is given in degrees; 
all other values must be given in units 
as specified by the flag IUNITS ' 

Flag determining the type of initial con- 
ditions input 

1 = Keplerian elements 3 

2 = Cartesian coordinates 

y. 

3 = Apogee and perigee 0 

Final condition that must be satisfied 
in order to stop program execution 
normally 

Flag that specifies the type of STOP 
condition ■ 

1 = STOP in days 

2 = STOP in revolutions 3 

Increment- for the printed output; 
a value is not needed if IPRINT is 
set to 0; PRINT = 0.0 is a valid 
entry 


a Default value. 

^To be implemented. 
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Variable 

IPRINT 

DATE 

IDRAG 

CD 

AREA 

XMASS 

ILONG 


a Default 


TABLE I.- Continued 


I 


DP 


No. of inputs 
required 


Description and 
available options 


1 


Flag that specifies the PRINT increment 


0 = No PRINT increment 3 

1 = PRINT is days 

2 = PRINT is revolutions 


6 


Date of epoch given as a calendar date 
of the form 


I 


1 


DP 1 

DP 1 

DP 1 


I 


1 


Month, day, year, hours, minutes, 
seconds 

Note: Computation range is from March 1, 

1900 through February 28, 2100 

Flag that specifies if the drag equations 
are to be included in the computation 

0 = No 

1 = Yes a 

Coefficient of drag, a value is not 
needed if IDRAG is set to 0 

Frontal surface area of the satellite; a 
value is not needed if IDRAG is set to 0 

Total mass of the satellite in kilograms; 
a value is not needed if IDRAG is 
set to 0 

Flag that specifies the type of poten- 
tial terms to be included in the computa- 
tions 

0 = None (two-body orbit) 

1 = J 2 short period, and first-order sec- 

ular terms 3 

2 = Compute the mean energy due to 

geopotential terms as defined by 
NMAX and MMAX 


value . 
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TABLE I.- Concluded 


Variable Type 
NMAX I 

MMAX I 

IPSPRT I 


IUNITS I 


No. of inputs Description and 

required ' available options 

1 Total number of zonal terms to be 

included by the geopotential model; a 
value is needed only if ILONG ’ is set 
to 2 

1 Total number of tes serai terms to be 

included by the geopotential model; a 
value is needed only if ILONG is set 
to 2 

1 Flag to determine if the PS elements are 

to be included with all printout 

0 = No a 

1 = Yes 

1 Flag that specifies the units of the 

input data and selects the appropriate 
physical constants 

1 = km, sec a 

2 = ran, sec 

3 = ft, sec 

4 = m, sec 

5 = km, hr 

6 = nm, hr 

7 = E.r., min 


a Default value. 
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2.1.2 Default Values 

To help shorten the number of data values that must be supplied by the user, 
the ASOP program assumes certain default values for those variables not 
explicitly mentioned on the input NAMELIST. These default values are listed 
in table II and a description of the variables can be found in section 2.1.1. 
Any variable not listed in table II must be specified by the user. 


TABLE II.- DEFAULT NAMELIST VALUES 


Variable 

Default value 

IEL 

1 

STOP 

100.0 

ISTOP 

2 

PRINT 

0.0 

IPRINT 

0 

DATE 

1 . , 1 . ,1978 

IDRAG 

1 

AREA 

185.3 

CD 

2.2 

XMASS 

90700.0 

ILONG 

1 

NMAX 

2 

.MM AX 

0 

IPSPRT 

0 

IUNITS 

1 
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2.1.3 Output Description 

After the ASOP program has been started with the command 
gXQT *NUMEG. ASOP-PROG 


and the input data 
NAMELIST variables 
in figure 2. 

$ INPUT 
EL = 

99D+004, 

99D+001, 

OOD+QOl 
I EL = 

STOP - 

I STOP - 

PRINT = 

I PRINT ~ 

DATE = 

0QD+Q01 , 

OOD+OOO, 

QOD+OQO 
IDPPG = • 

CD = 

APEfi = 

::mrss = 

I LONG = 

MM A X 

mma;-; 


has been added, the program will print out all of the 
as shown in figure 1 and the initial conditions shown 


. 3SH207O5H999999999D+S0H. 

. 2152207510000000000+004-, 

. 31 lSH-50S0000000O00D+00i , 
4- 2. 

. -100OO00000000O0000D + 003 
+ 2 

_ 5OO0GOOSOSOQQOOO0CID + QO2 
+ 2 

. 1000000000000000000 + 001 , 

. 1978000000000000000+004-. 

. 0000000000000000000+000. 
+ 1 

. 2200000000000000000+001 

. 1853000000000000000+003 

. 9070000000000000000+005 
+2 
+4 
+ 4 


IPSPPT = 

+ 1 

I UN ITS = 

+ i • 

TEND 


Figure 1 NAMELIST data output 


INITIAL CONDITION 1 : 

Cl/01/1978 00:00: 

. 000 

DAYS: 0.O00U000 

! REVS: 0.0000000 

fl= 5.6993532+03 

KM ! E= .0009100 

0MEGP= 17.92236 

DEG !N0DE= 20-OQO00 


5255 1 78579999999 
54 1 478 3099999999 
2955928820000000 


, 1000O0O000O0000O 
0000000000008000 
0000000000000000 


2443509.5000000 
• + (CHECK: 0.0000000 
— 

30.00000 DEG 

22.05838 DEG 

2. 1522075 + 03 i M . 
2.9559288+OQ IM.-'S 


v::= 

PS ELE 


-5.4147831+00 KM/Si 
rtENTS: 

.10475457+O01 . 59495155-0O2 

.51675514+005 .86357683-001 


! 1 = 

, DEG ! M= 

i— 5. 2561786+03 KM ! Z— 
VY= 3.1154505+00 KM/S!VZ= 


>, 

!• I 


-.40245815+001 
. 11057447+O03 


ENTER: :{= EXECUTE; S= STOP: 


' '■ : I : *: ! J* : :Js •• { 

C= CHANGE INPUT 


9373237+0OO 
.29759353+002 

i ! ! ! 


Figure 2.- Initial condition format. 


After the initial conditions have been displayed, the program will wait for 
the user to check the input parameters. Some cheeks that can be made are 
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a. CHECK value should be 0.0000000. 

b. DAYS and REVS should be 0.0000000. 

c. A double asterisk (**) will appear after the output condition that is 
being satisfied, i.e., after the DAYS or the REVS value. This will, agree 
with input value of IPRINT or IST0P. 

d. Are the initial conditions (a, e, i, etc.) the desired values? 

e. Are the units correct (km, km/s,etc.)? 

Once the user is satisfied that the printed initial conditions are the ones 
needed, an X should be entered. When the X is entered, the program will 
continue execution and print out information as specified by the input parameters, 
i.e., STOP, IST0P, PRINT, AND IPRINT (figure 3). 


ENlEh-': .;= EXECUTE; S= STOP; C= CHANGE INPUT 


A SOP OUTPUT 

01/04/1978 03=38=59.070 

DA VS : 3.1520726+00 ! REVS: 


0000000+01 


FREE: 0002 
2443512.6520726 
** ! CHECK : -3.7947667-06 


R- 6.6965691+03 

H 

KM 

h 

! E= 

-00US5H 

! 1 = 

29. 9362S 

DEG 

0MEER= 31.71316 

DEG 

!N0DE= 

366.98032 

DEG ! M= 

30.94208 

DEG 

:;= 3. 332954 4+0 3 

KM 

i Y= 

4 . 98 1 H 775+0 3 

KM ! Z= 

2. 9717888+03 

KM 

v: 6907 30 3+00 

KM/S ! 

! VY= 

3 . 4-262523 + 00 

KM/S ! VZ= 

1 . 770928O+-O0 

KM/S 


. 307411 38+002 
.51664876+005 


47544951-001 
1344 3.374+000 

i V V 4 Jv{ | \ :j 


-61953*16+001 
. 117^277+003 


2723^011+006 
297735 3 1+002 


K 


01/07/1973 07=17'; 
DAYS: 6.3034791+00 


FINAL CONDITIONS 

5^6 

! REVS: 1.0003000+02 «= H CHEC 

+ + 

fi= 6.6946143+03 KM ! E= .0014630 i I=- 

OMEGAr 48.11705 DEG ! N0DE= 333.^3645 DEE i M= 

::= 3 . oi 36710+03 km \ v= 4.9535461+03 1 m i z= 

V:;=-6. 6812825+00 KM/S! VY= 3.S65727O+O0 KM/S!V2= 

PS ELEMENTS: 

.32717351+002 -.76850190-001 
. 51657340+005 . 27972934+000 

i : :! #e . ;•:! s{: :!:•!: } l :!•*:} 

INPUT DATA USING THE NAMELIST ' TINPUT-' 


244351 

: 9.Q 


5.30347 

538672- 


91 

08 


29.97883 DEG 

37.22867 DEG 

3.3307289+03 KM 
3.0931351-01 KM/S 


.51655997+00 
. 1O561303+OO 
* : *! 


3 


m SHHS2 
- 29777 
* 1**1 ['• 


213+906 
608 + 002 


Figure 3*- Intermediate and final output format. 

ASOP will always print the conditions as stated by PRINT and STOP unless it is 
unable to do so because of satellite reentry. The reentry condition is defined 
as the point where ASOP determines the satellite will reenter within one revolu- 
tion. A sample of reentry is shown in figure 4. 
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12/02/1978 00:00: -000 

DAYS: 0.0000000 ' ■ 


INITIAL CONDITIONS 
! REVS: 0.0000000 


2443844. 50OO00O 
i CHECK: 0.0000000 


f|— F. ^201400+02 KM ! E= .0000000 1 = 

GMEGA= .00000 DEE ! NODE = uO-QOOOO DEG M- 
g, 1 ^8° 274+03 KM ! 7= 2--2.3QO-192+-0.J -} M ! 

Vil=-2;34789S'O+0O KM/S! VY= 6. 4507885+00 KM/S' VZ= 

pc, fI'fmEWTS: 

xCo0 t c 85+rjO o .79599283-001 -.38148 354+002 - 
! 5097979 i +005 -.21889718+000 .10481x74+003^ 

ENTER: EXECUTE; S= STOP; C= CHANGE INPUT 

“aSOP OUTPUT 


28.80000 DEG 

.00000 DEG 

0 .0000000 KM 

3.7428027+00 V M/S 

.80044931-018 
. 30598443+002 
: : ; s+.: r-;= - : : 


PAGE: 0002 


12/02,- 

DAYS: 


RENTPY CONDITIONS 

- - 7 CT - £ ^ 7^0 

507-82 " " 1 PEVS= 1.0946538+00 


-MCHEC 

* i I = 
i M = 


244 3844. 566154-5 
K: 2- 9 364986-06 


OKEG fi- 


PS tLE 


3 - 9 t 
-6- OS 
MENT« 


508*4 

1.C72 


14+03 PM 
66 DEG 
2S+03 KM I 
00+00 KM/S! 


i E= 
I NODC= 


. 0007424 
19.49643 DEG 
4 . 8441 56 3+0 3 KM 
3. 8439442+00 KM/ 


^ ' T 

INPU'’ 'i 


: 4 *«* 337 3 + 000 -.17842732-001 - 
1:0934787+005 -.13947939+000 

i ; ,V'| * • :*- J ■ J. i , ’! t :!*:!: ! 

■i USING THE NAMELIST ’flNPUT 1 


37192987 

10505058 


+002 

+003 


28.58721 DEO 

312.57730 DEE 

1.7870120+03 KM 
3.0831513+00 KM/S 

.57148404+004 
. 30845447+002 


Figure 4.- Reentry format. 


In some cases, ASOP may determine that the satellite will reenter within 
one revolution of the initial conditions. Therefore, ASOP would print out 
the initial conditions as the reentry conditions. 

In general, all output is clearly labeled, but some terms should be explained 
further . 

DAYS: Total number of days elapsed since the starting epoch. 

REVS: Total number of revolutions completed. 

CHECK: Value indicating the accuracy of the analytical theory; although 

this value is necessary as a check on the theory, it is not 
a sufficient check. 

PS ELEMENTS: The Poincare-Similar elements listed as 

a 1 °2 a 3 a 4 
Pf P2 P3 P4 

Double asterisk (**): Indicates, the stopping condition being satisfied*; 

> this flag will move between the DAYS and REVS 
value as needed* 


12 



77FM50 


2.1.4 Run Setup (Control Cards) 

The ASOP program is written in standard FORTRAN V and designed to run on 
the NASA/JSC UNIVAC 1110 computer using the EXEC-8 operating system. All 
the relocatable and executable elements are on the file FM6-N08569*NUMEG. 
ASOP may be executed by entering the following for demand operation. 

a. gQUAL FM6-N08569 

b. gASG.A *NUMEG. 

c. gXQT *NUMEG. AS0P-PR0G 

d. Add input data 

e. Enter the letter X, S, or C (see section 2.1.3) 

f. gEOF or go to step d. 

g. gFIN 

If run in a batch mode, the following input cards are needed. 

a. gQUAL FM6-N08569 

b. gASG, A *NUMEG. 

c. gXQT *NUMEG . AS0P-PR0G 

d. Add data file or data cards 

e. Card containing the letter X in column 1. 

f. Repeat instructions d and e as often as necessary. 

g. gEOF 

h. gFIN 
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2.1.5 Sample Computer Run 

In this, section, a sample computer run is reproduced for a typical Shuttle-type 
orbit. The orbit has been predicted for 100 revolutions ( = 6.3 days) with the 
output given every three days. 

This example is intended to familiarize the user with the format of the ASOP 
output and to illustrate the use of the various input options discussed in sec- 
tion 2.1.1. A full description of the output is given in section 2.1.3. - 

Initial parameters 


Semimajor axis (a) 
Eccentricity (e) 
Inclination (i) 

Argument of perigee (w) 
Argument of the' ascending 
Mean anomaly (M) 


6699.3532 km (1.05 ER) 

.001 

30 degrees 
18 degrees 
(ft ) 20 degrees 

22 degrees 


Input parameters under the NAMELIST $INPUT 
$INPUT 

EL(1) = 6599.3532 

EL(2) = .001 

EL(3) = 30.0 

EL(4) = 18.0 

EL(5) = 20.0 

EL(6) = 22.0 

IEL = 1 

STOP = 100.0 

IST0P = 2 

PRINT =3.0 

IPRINT = 1 

DATE = (Default Values Used) 

IDRAG = 1 

CD = 2.2 

AREA = 185.3 

XMASS = 90700.0 

ILONG = 2 

NMAX = 8 

MMAX = 8 

IPSPRT = 1 

IUNITS = 1 

$END 


14 



77FM50 


Sample computer run. 


. KNOT -HUMES- ASOP-PEOG 
INPUT DATA USING THE NAMELIST 'f INPUT' 
? KfiDD *NUMEG . DATA-RS0P/L0NG2 


f END 




* INPUT 

EL 

- 

. 6699 353 19999999998D+004, 

_ 9999999999999999 

97D-003, 



- .30OOO0OOOOOOO0OO0OD+002. 

- 1800000000000000 

OOD+OGll, 



-2GOOOO0OO00OOOO0O0D+002, 

. 2200000000000000 

00D*C02 

I EL 

r 

+ 1 


STOP 

- 

. 10QOGG0O08GGOG00GQD+OO3 


ISTEF 

~ 

+2 

„ 

PRINT 

JZ 

. 30000O000O00O000OOD+OO1 


I PRINT 

ZL 

+ 1 


DHTE 

z 

. 1000000000000000000 + 001 . 

_ 1090000000000000 

O0D + UO 1 , 

GGD+008, 


. 1978OOO0OOOOOOO0OOD+004, 
. OOO0Q0QOOQO00OOO0OD+00O, 

.0000000000000000 

.0000000000000000 

OCD+OO0 

IOPRL 

“ 

+ 1 


CD 


. 22000OO0OOOO00O8OGD+001 


RPER 


. 1 S5300GG00GG0OGO0QD+00 3 


i 1MRSS 

— 

. 9O7GG0OGGOGGG0GO00D+GGS 


IL0N13 

z 

+2 


who:: 

“ 



mmm: ; 

z 

+ 8 


IPSPRT 

z 

+ i 


I UNITS 

z 

+ i 


*END 


01/01/1970 00 : 00 : .000 

DAYS: 0.0000000 


INITIAL CONDITIONS 
! REVS: O.OO0OO0O 


2443509. 5QQOO0O 
'CHECK: 0.OO0OOOO 


S = 6.6993532+03 KM ! E= .0010000 ! 1= 

OMEGA= 18.OOO0O DEG ! NODE= 20.OOOOO DEG ! M= 

2= 3 5 395374+03 b M ! Y= 5.2568222+03 KM I 7= 


VY= 3. 1134747+00 KM/S ! VZ= 2.95626Q8+0Q h M/$ 


V;;=- 6 . fcd287+00 HM/S! 

PS ELEMENTS: 

.10479476+001 -.76954051-002 -.40245812+602 
.51675611+005 .16285596+000 .11057446+003 


30.00000 DEG 

22.00000 DEG 
2. 1530569+0 3 b M 


. 65091474+080 
. 29769235+082 


-I * <s !- ! 1**1 I -I *!- i l I* *5* *1 ) ! =!* I I ! ! •! •!- ”! ! *1 : ! : ! I *1 •!• « : I I !" I I ■!* “!* ! ^ S' !* * ! •!'•!■*}* I * 'I *"'*!* !* I ! !••!• !• I { i '» 

ENTER: ;;= EXECUTE; S= STOP; C= CHANGE INPUT 
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Sample computer run.- Concluded 


ENTER: X= E/EGUTE ; S= STOP; C= CHRNGE INPUT 

B I 

"fiSOP OUTPUT P'flGE: 0002 

01/04/1978 00:00- .000 24H3512.500OOOG 

DAVj : 3.0000000+00 + : I PE VS: 4.7586650 + 01 ! CHECK: 1.4409285-05 

J 1 

A = 6.6955019+03 J- M ! E= .0010075 ' 1= 29.97867 DEG 

OMEGA:: 318. +■3280 DEG iNODEs -358.0-5518 DEG i H= 314.47288 DEG 

.)= 1. 55^4173+02 KM ! V=-5. 7964700+03 KM ! Z=-3. 3390751+03 kM 

V::= 7.7162475 + 00 KM/S> VV= 7.8156721-02 1- M/S 1 VZ= 1.9615716-01 1 M/8 
PS ELEMENTS. 

.1-357814 7+002 .10185805+000 . 39900710 + 001 .25919875 + 006 

.51661116+005 . 16527396+000 .11750491+083 .29773243+802 

t t t ( i > 1 . t i r r .1 -i. . . ■ r • • i t i .1 r i i .1. I t ».«•■.» 1 t t i« I* 1 > t .1 .1. 1 I. .1 t . 1 . r 1 1 1 1 1 1 1 1 r I 1 » I .1 .1 » 1 .« 

14 I « I I I I I »■*♦**! I 1 I I t I t * < II I *1 II I I I • • 1 *1 ' I' I' 1 !• 1 ■ • ■ ■ I I I* ' » I I *1 «' I I I I* l I lit! ■ I I I •(• 

INITIALIZATION TIME = 569 MS. 

01/07/1978 00:00: .000 2443515.5000000 

D H V S : 6. 0000000 + 00 =!=:- ! REVS: 9.51339-35+01 ! CHECK: 1.1420590-85 

+ + ^ 

A= 6. 698-3353-^0 3 KM ! E= .0011075 ■ 1= 38.00575 DEG 

OMEGAS 88.20042 DEG * NGDE = 3-36.11077 DEG ! M= 61.22178 DEG 

:;=-4. 085^234+03 KM ! V= -5.0245406 + 0 3 KM 1 Z= 1.6975347+03 KM 

VK=-5. 9157146+00 KM/S! VV=-3 . 676952 3+00 h M/S ■ VZ=- 3. 3249412+00 KM/S 
PS ELEMENTS: 

.52458389+002 82058586-001 .47660244+002 .51840169+006 

.51671700+005 .21288088+000 .10760608+003 .29777083+002 

* J -I* fl l f 1 T T» 1* «■«. Til «l« 1 I •!« *1 T I «I I* I t t t 1**1* I T «f« * I * » I * l« I I t » T I* »l I *1* P »f T -I | f f 1**1 •!* •■»*!* 1 T I. K f 

« i* A 4 *■* r*r « i i a a i ••• a i i*n« A i i**r i i* < a a 4 M 1 A a 'a* *1 r r » » i- i *p t i i «j i 4 i *1 i* 4 i 1' 4 


FINAL CONDITIONS 

01/07/1978 07:17: 7.391 2443515.8835578 

DATS: 6. 3035518+08 ! PEVS: 1.0800000+02 : !”i > CHECK: -5. 9225315-06 

“ ■“ ~~~ 1 ! * — — — — . — — 

A= 6.6946143+03 I'M ■ E= .0014899 i 1= 29.97886 DEG 

OHEGA= 51.13169 DEG ! NODE= 333.93334 DEG ! M= 34.24363 DEG 

3.011288 3+03 KM ! V= 4.9544626 + 03 IM ■ Z= 3. 3306205 + 03 IM 

V:i = -6. 6831533 + 00 K M/S ! VV= 3.8636606+00 V M/S l VZ= 3.0302620-01 f M/S 
PS ELEMENTS: 

-.18745875+003 - . 93856906-00 1 .51661785 + 002 .54462885+806 

. 5a657337 + 0C5 .27737-512 + 000 .10561032+003 . 29777855 + 002 


INPLT DATA USING THE NAMELIST 1 -f INPUT’ 
%EGF 

: : NORMAL PROGRAM TERMINATION : : 
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In the interactive mode, the program solicits input data with the messages: 

a. INPUT DATA USING THE NAMELIST '$INPUT* 

b. ENTER: X = EXECUTE; S = STOP; C = CHANGE DATA 

c. INPUT NEW DATA USING THE NAMELIST '$INPUT' 

In the batch mode, the data input must be followed by one card containing the 
letter X as shown in the following examples: 

Example (1 ) Example (2) 

b$INPUT gADD filename -data element 


necessary data cards changes to data element 

(see section 2.1.1) (if any) 


K$END 

X 

b$INPUT 


tt$END 

X 

§ADD filename . data element 


etc. 


etc. 


• • 

©EOF @EOF 

@FIN @FIN 

(b is a blank that must be included) 


2.2 SUBROUTINE PACKAGE 

Along with the ASOP stand-alone program, there is a subroutine package that may 
be included in the user's software. This package is in the form of a relocat- 
able element and is located in FM6-N08569*NUMEG. ASOP-SUB. Its operation is 
identical with any user-written subroutine. 

This section will describe the information needed by the user to ensure proper 
insertion and operation of the ASOP package within the user's software. 
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2.2.1 Required Subroutines 

The ASOP subroutine package consists of 27 subroutines. • These are a driver 
subroutine (ASOP) that controls the basic logic of the package, 10 general 
subroutines that perform the functions necessary to .the analytical theory 
and 17 subroutines to. initialize the drag model and the computation 
of mean energy. 

These subroutines are 


ASOP 
CONST 
DRAG 
GEOPOT 
LON GPP 

POTEXP 

PREPD 

PSANS 

PSTOX 

TIMEPS 

XTOPS 


Driver subroutine 

Planetary and mathematical constants 

Adjust the PS elements to account for drag perturbations 
Determine Earth's gravitational potential 
Compute first-order zonal long periodic perturbations and 
second order zonal secular perturbations 

Compute mean energy due to tesseral and sectorial geopotential 
harmonics 

Initialize the drag model 
PS analytical Jg theory 

Transformation subroutine: PS elements to Cartesian coordinates 

Time iteration stopping procedure 

Transformation subroutine: Cartesian coordinates to PS 

elements 


Initializing subroutine called by CONST: 

PREPT Initialize the geopotential coefficients for the Earth 


Computational subroutines called by LONGPP 


DETERM Compute the first-order long-period generating function 
and its derivatives 

PPRIME Compute the second-order zonal Hamiltonian and its derivatives 

Initializing subroutines called by POTEXP: 

COEFF Compute the binomial coefficients and the Fourier coefficients 
of the powers of cosine and sine 

IL0G10 3 Determine the number of terms to be included in the expansion 
for the Earth’s geopotential model 

RECUR Compute the sine and cosine of multiples of an angle 

TIMEXP Compute coefficients of the expansion of the time equation 


a IL0G10 is also called by DETERM and FPRIME. 
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Initializing subroutines called by P8EP: 

Compute the canonical forces due to atmospheric drag for 
the PS equation 

Compute the atmospheric density at a given altitude above 
an oblate Earth 

Initialize the coefficients for the Jacchia 71/Lineberry 
atmospheric density model 
Invert an n x n matrix and/or solve Ax = B 
Convert .the mean anomaly to eccentric anomaly and compute 
its sine and cosine 

Establish the parameters needed to calculate- the position 
of the Sun 

Determine the solar activity coefficients for a given date 
Compute the position of the Sun analytically 
Generate the table of coefficients for the sine and cosine 
of tfJ i for the atmospheric drag function 

The subroutines listed above are fully described in section 3.3, and a diagram 
of the data flow between these subroutines can be found in appendix H. 

To help the user add these subroutines to his own software, a relocatable ele- 
ment has been formed that includes all the above subroutines. Therefore, the 
user needs only to include the element 


FM6-N08569 *NUMEG . ASOP-SOB 


when forming an executable element. 


2.2.2 Interface Requirements 

To access the ASOP subroutine package, the programer must use the FORTRAN state- 
ment 


CALL ASOP ( X , STOP , ISTOP , NEWX ) 


A full description of the argument list variables can be found in sections 2.2.3 
and 3-3.3. Also, the user must initialize certain COMMON block variables before 
entering the ASOP subroutine. 

Table III gives a list of the variables that must be initialized prior to call- 
ing the ASOP subroutine and the COMMON block in which the variables are located.. 
A complete description of the variables and their allowed values can be found in 
section 2.1.1. 


CANFOR 

DENSTY 

INITAL 

MATIN 

MTOECC 

PREPS. 

SACT 

SUN 

TABLE 
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TABLE III.- 

COMMON BLOCK 

INITIALIZATION 

Variable 

COMMON 

block 

Default 
value (s) 

IUNITS 

CPRINT 

1 

IDRAG 

PERTRB 

1 

CD 

DRAG 

2.2 

AREA 

DRAG 

185.3 

XMASS 

DRAG 

90,700.0 

DATE 

EPOCH 

•1.0,1.0,1978.0 
0.0, 0.0, 0.0 

XJDATE 

EPOCH 

2,443,509.5 

ILONG 

PERTRB- 

1 

NMAX 

TESS 

2 

MMAX 

TESS 

0 


2.2.3 Input/Output Description 

The argument list to the ASOP subroutine consists of four arguments given in the 
following order. 


CALL ASOP (X , STOP , ISTOP , NEWX ) 


On input the variables are 


X An array of eight elements corresponding to the initial state vector 

in the following order: 

X(1): X position component 

X(2): Y position component 

X(3): Z position component 

X(4): X velocity component 

X(5): y velocity component 

X(6): Z velocity component 

X(7): Physical time (set to zero) 

X(8): Total energy (set to zero) 

STOP Stop value desired; it may be given in days or revolutions. 
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WARNING: A number must not be used in the following positions of the argument 

list. Assign the desired value to a variable, and use the variable 
in the argument list. If a number is used in these positions instead 
of a variable, unpredictable results may occur. 


ISTOP Flag determining whether the value given to STOP is in days (ISTOP 
= 1) or revolutions (ISTOP =2). 


WARNING: Only a 1 or 2 should be used as input. On output, ISTOP should be 

checked to see if ISTOP was reset to 3 indicating reentry conditions. 


NEWX Flag determining if the ASOP subroutine is to be initialized 
NO = 0, YES = 1 

The initialization process must be done whenever new initial condi- 
tions are entered. 

Input to the ASOP program is also done by means of COMMON blocks. These COMMON 
block input variables control the internal operation of the ASOP subroutine pack- 
age and should not be changed once the subroutine has been initialized. 

Table III (section 2.2.2) gives a complete list of the variables that must be 
initialized and their default values. 

On output the variables are 

X An array of eight elements corresponding to the final state vector at 
the given value of STOP. The order is the same as for input. If the 
value of STOP was given in days (ISTOP = 1), then the value of X(7) will 
be an approximate value of STOP given in the time units specified by 
IUNITS . 

STOP Unchanged from the input value. 

ISTOP Set to three (3) if reentry condition exists. Otherwise, it is 

unchanged from the input value. Therefore, a variable name should 
always occupy this position in the argument list. 

NEWX Set to zero (0) upon return from the ASOP subroutine; therefore, a 
variable name should always occupy this position in the argument 
list. Otherwise, unpredictable results may occur. 
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3.0 DESCRIPTION AND STRUCTURE OF ASOP 


The ASOP program was designed as a top-down structured program consisting of 33 
subroutines (modules) and a main (driver) program. Within this set of subrou- 
tines, a subset of 27 subroutines comprises the removable ASOP subroutine pack- 
age . The. package has been designed for easy incorporation into existing user 
software . 

The two sets of subroutines are described in the following subsections and each 
subroutine is documented in section 3-3. 


3.1 THE ASOP PROGRAM 

Basically, the ASOP program, shown in figure 5, consists of four segments: a 

main program or driver , an input routine , an output routine , and a removable 
ASOP subroutine package. The removable package is described in section 3-2 so 
that only the first three segments will be described here. 

The purpose of the main program is to call all the necessary subroutines (input, 
output, constants, etc.) in a specific sequence to produce the desired results. 

In particular, the main program provides for the repetitive call to the ASOP 
subroutine in order to produce a satellite ephemeris . 

All input to the program is controlled by the subroutine INPUT. Its primary 
functions are 

a. Set all default values 

b. Accept input from the NAMELIST $INPUT 

c. Convert any input values that require conversion 

d. Issue normal program termination command 

Output from the ASOP Program is performed only by the subroutine OUTPUT. This 
subroutine contains all the FORMAT specifications used for normal 3 output from 
ASOP. The output program will also perform any conversions required to make the • 
output more understandable. This involves converting Cartesian coordinates to 
Keplerian elements, radians to degrees, and time units to days. 


a All error output is controlled by the individual subroutines. 
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Page 1 of 1 

Figure 5.- Structure of the ASOP program. 
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3.2 THE ASOP SUBROUTINE PACKAGE 

The ASOP subroutine package has been designed as an independent segment of the 
ASOP program so that it can be incorporated into existing software with little 
or no modifications. 

This package has seven basic parts (fig. 6): 

a. A driver subroutine (ASOP) 

b. Coordinate transformations (XTOPS, PSTOX) 

c. A stopping routine (TIMEPS) 

d. The analytical theory (PSANS, LONGPP, DRAG) 

e. Initialization of the geopotential model 

f. Initialization of the drag model (PREPD) 

g. Initialization of physical and trigonometric constants (CONST, PREPT) 

Along with the 11 subroutines mentioned above there are an additional 16 subrou- 
tines. Subroutine GEOPOT performs computations required by the coordinate trans- 
formation subroutines XTOPS and PSTOX. Subroutines LONGPP, POTEXP, and PREPD 
use the remaining 15 subroutines for initialization of their respective 
perturbation models. 

Subroutine ASOP performs the same task as the main program in that it calls all 
the necessary subroutines as dictated by the input values. In the case of the 
ASOP subroutine, however, the input values are given through an argument list 
and a few COMMON blocks. Therefore, if a satellite ephemeris is desired, the 
user must supply the necessary coding within his own software. 

The basic input to the subroutine is: 

a. The Cartesian coordinates (X) 

b. A stop value (STOP) 

c. A stop flag (ISTOP) 

d. A new data flag (NEWX) 

This input is fully described in sections 2.2.3 and 3.3.3 

Because the analytical theory has been developed in PS elements, it is necessary 
to perform transformations to and from the element set. Transformation into the 
PS elements from Cartesian coordinates is performed by the subroutine XTOPS 
while the reverse transformation is performed by the subroutine PSTOX. 
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A time-stop subroutine has been included because the PS elements use the true 
anomaly as the independent variable and do not use time. This stopping routine 
is an iterative procedure and is described in section A. 3. 

Subroutine PSANS updates the PS elements at a specified value of the independent 
variable. The J 2 perturbations are computed within PSANS, the drag pertur- 
bations are computed by a call to subroutine DRAG, and the long-period geo- 
potential terms are computed by LONGPP. 
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Figure 6.- Structure of the ASOP- subroutines. 


26 














77FM50 


3-3 MODULE DESCRIPTIONS 

'piis section will give a complete description of the subroutines currently used 
in the ASOP program. Each description will contain a brief statement as to the 
purpose, or use, of the subroutine as well as a description of important vari- 
ables used within the subroutine. Also included are lists of the named COMMON 
blocks used, external references to other ASOP subroutines, and other ASOP sub- 
routines that reference the subroutine being- described. Information is also 
available as to the calling sequence of the subroutine and the size of 'the 
subroutine. Each description is followed by a general flow chart of the sub- 
routine (figs. 7 through 40). 

Each program is listed alphabetically, with the exception of the MAIN program, 
which is described first. 


3*3*1 MAIN P r ogram (Driver) 

Purpose : Driver for the Analytical Satellite Orbit Predictor (ASOP) program 

Calling sequ ence: None 

Called by : Operating system 


Subroutines /functions called: 


Named COMMO N: /CARTC / 

/ CBASIC/ 
/CPRINT/ 
/END / 
/PS / 
/PSTIME/ 


ASOP, INPUT, OUTPUT 
X(8) ,R,RI 

PI , TWOPI , DEG , RAD , D AY , DTOKM 
PRINT , IPRINT , IPSPRT, IUNITS 
STOP, ISTOP 

SIG{ 8 ) , TAU , TAUMAX , TAUINT 
CLO , FAKTS , TOL 


Equivalence : (X(7),TIME) 

Program dat a: Size = 1478 C 1 03 1 o ) wor> ds compiled. 


FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

DAY 

1 

DP 

I 

Conversion of days into hours, 
minutes, or seconds. Its numeric 
value is controlled by IUNITS 

IPRINT 

1 

I 

I/O 

Flag to determine if the inter- 
mediate printout is to be given 
on days or revolutions. If re- 
entry condition exists, flag 
is reset to 3. 
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FORTRAN Input/ 


variable 

Dimension 

Type. 

output 

Description 


• 



= 0 no printout 
•= 1 days_ 

= 2 revolutions 
= 3 reentry condition exists 
(output only) 

ISTOP 

1 

I 

I/O 

Flag to determine if the 
program is to terminate 
after a given number of 
days or revolutions as 
specified by STOP. If 
reentry condition exists 
flag is reset to 3 





= 1 days 
= 2 revolutions 
= 3 reentry condition exists 
(output only) 

NEWX 

1 

I 

I/O 

Flag to determine if the ASOP 
program is to be initialized 





= 0 no 
= 1 yes 

PRINT 

1 

DP 

I 

Increment for which the inter- 
mediate printout is desired 
(valid only if IPRINT >1) 

STOP 

1 

DP 

I 

Value at which the program is 
to stop execution (units are 
determined by ISTOP) 

STOPDT 

1 

DP 

0 

Value at which the next 
intermediate printout is 
desired (valid only if 
IPRINT > 1 and PRINT > 0.0) . 

TAU 

1 

DP 

I 

Independent variable of 
the PS elements; it is 
defined such that REVS 
=• TAU/2TT 

TAUINT 

1 

DP 

I 

Initial value of TAU for which 


the force model is valid 
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FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

TOL 

1 

DP 

I 

Tolerance used by the 
time- stopping routine. 

TWOPI 

1 

DP 

I 

2TT- 

X 

8 

DP 

1 

Space vehicle parameters 


in the form: 

XC 1) + X(3) = X 
X(4) -«■ X(6) = V 
X(7) = time 
X(8) = total energy 
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3-3-2 AEIXYZ (Subroutine) 


Purpose : Transform the Keplerian elements (a,e,i,oi,fl,M) into Cartesian coordi- 

nates Of,V) 


Calling sequence : CALL' AEIXYZ 

Called by: INPUT 


Subroutines/ functions used : MTOECC 


Named COMMON: 


/CARTC / 
/CBODY / 
/KEPLER/ 
/RPOOL / 


X(3),V(3), TIME, ENERGY, R,RI 

XMU , XMUI , SQTMU , SQTMUI , EPS 

A , E , XI , OMEGA , XNODE , XM 

SINC,SOMEGA ,SNODE , CINC , COMEGA , CNODE , 

B1(3),B2(3),X11,X12,V11 ,V12, TEMPO, TEMPI 


Equivalence : (XI ,XX0( 1 ) ) , (SINC,XX1 ( 1 ) ) , (CINC',XX2( 1 ) ) 


Program data : Size = 302$ (194-) 0 ) words compiled 

Subroutine valid only when | e| <1.0 


FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

A 

1 

DP 

I 

Semimajor axis of the 
satellite’s orbit; must be 
greater than 0 

COSEA 

1 

DP 

I 

Cosine of the eccentric 
anomaly (cos E) 

E 

1 

DP 

I 

Orbital eccentricity (e); 
must not be greater than 1 

EA 

1 

DP 

I 

The eccentric anomaly of the 
satellite computed from 
Kepler’s equation (rad) 

OMEGA 

1 

DP 

I 

Argument of pericenter (o) 
(rad) 

R 

1 

DP 

0 

Magnitude of the position 
vector of the satellite 

RI 

1 

DP 

0 

Inverse of R 

SINEA 

. 1 

DP 

I 

Sine of the eccentric anomaly 
(sin E) 
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FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

V 

3 

DP 

0 

Velocity vector of the satel- 
lite with respect to the 
Earth's equator 

V(1) = v x 
V(2) = Vy 
V(3) = T z 

X 

3 

DP 

0 

Position vector of the satel- 
lite with respect to the 
Earth's equator 

X(1) = X 
X(2) = X 
X(3) = Z 

XI 

1 

DP 

I 

Orbital inclination to the 
Earth's equator (i,rad) 

XM 

1 

DP 

I 

Mean anomaly of the satellite 
(M,rad) 

XMU 

1 

DP 

I 

Gravitational constant for 
the central body (U) 

XNODE 

1 

DP 

I 

Argument of the ascending 
node (&,rad) 
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Figure 8,-AEIXYZ flow chart. 
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3.3.3 ASOP (Subroutine/Driver) 


Purpose : Driver for the analytical section of the ASOP program; it controls all 

the operations required by the analytical program PSANS 

Cal-ling sequence : -GALL ASOP (X , STOP , ISTOP , NEWX ) 


Called by : MAIN 


Sub routines/ functions used : CONST 3 , POTEXP, PREPD’, PSANS, PSTOX, TIMEPS, XTOPS 


Named COMMO N: 


/CARTC / 
/CBASIC/ 
/GEO / 
/PERTRB/ 
/PS / 


XIN(8),R,RI 

PI , TWOPI ,DEG , R AD , DA Y , DTOKM 
RE 

IDRAG.ILONG 

SIG( 8 ) , TAU, TAUMAX , TAUINT 


Program data : Size = 256q (174 10 ) words compiled 


FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

DAY 

1 

DP 

I 

Conversion of days into hours, 
minutes, or seconds 

DTOKM 

1 

DP 

I 

Converts distance into kilometers 

IDRAG 

1 

I 

I 

Flag to determine if the drag 
calculations are to be included 

= 0 no 

= 1 yes 

ILONG 

1 

I 

I 

Flag to determine the type of 


geopotential terms to be included 



= 0 

none (two-body orbit) 


= 1 

J 2 short period, and first-order 
secular terms 


= 2 

Compute the mean energy due to 
geopotential terms as defined 
by NMAX and MMAX (see table I for 
a description of NMAX and MMAX) 

a Called only in subroutine package. 
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FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

. Description 

ISET 

i 

I . 

I 

Flag to determine if the force model 
must be updated 





= 0 no (Force model still valid)' 

= 1 yes (Update state vector to 

T max an< * reinitialize force 
model) 

ISTOP 

i 

DP 

i/o' 

Flag to determine if the value of 
STOP is given as days or revolutions; 
if reentry conditions exist, flag 
is reset to 3 





= 1 days 

- 2 revolutions 

= 3 reentry 

NEWX 

i 

I 

I/O 

Flag to determine if the ASOP program 
is to be initialized 





= 0 no 

= 1 yes 

R 

i 

DP 

I 

Magnitude of position vector 

RE 

i 

DP 

I 

Equatorial radius of the central body 
(Earth) 

RI 

i 

DP 

I 

Inverse magnitude of the position vector 

STOP 

i 

DP 

I 

Value at which the final state vector 
is required; units are set by ISTOP 

TAU 

i 

DP 

0 

Independent variable of the PS elements 

TAUINT 

i 

DP 

I/O 

Initial value of t for which the 
force model is valid; initially set 
to 0 

TAUMAX 

i 

DP 

I/O 

Range of validity for the force model; 
when t exceeds this value ( T max )> 
the force 'model must be reinitialized 

TWOPI 

i 

DP 

I 

2TT 
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fortran 

variable 

Dimension 

lyp.e 

Input/ 

output 

X 

8 

DP 

I/O 

XIN 

8 

DP 

I/O 


Description 

Initial/final state vector 

X(1) + X(3) = X 
X(4) + X(6) = V 
X(7) = time 
X(8) r total energy 

If initializing (NEWX =1), then 
X(7) and X(8) will be set to 0 

Identical to X but allows the 
ASOP subroutine to be removed 
from the stand alone program 
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Figure 9.- ASOP flow chart. 
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Figure 9.- Concluded 






. 3. 3.4 CANFOR (Subroutine ) 


77FM50 


Purpose: Compute the canonical forces due to atmospheric drag for the PS equa- 

tion used with the ASOP program (refs. 6 and 7) 

Calling sequence : CALL CANFOR (CFORCE) 

Called by : PREPD 


Subroutines/ functions used: None 


N amed COMM ON: /CBODY / 

/DBETAS/ 
/ DTABLE/ 
/ GMTROT/ 
/PS / 
/PSANS1/ 
/PSANS2/ 
/RETRO / 


XMU , XMUI , SQTMU , SQTMUI , EPS 

B,BSQ,B3,B4 

T( 12) 

WE.THETAO 

.S1,S2,S3,S4,R1,R2,R3,XL,TAU 

SSIG 1(2), TWOL , XIQL , FAK 

SUM2 , SUM3 , BIFF 2 , DIFF3 , GC , HC , PSSQRT , PS , QS 

IRO 


Program data : Size = 524g (340 -|q) words compiled 


FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

B 

1 

DP 

I 

BSQ 

1 

DP 

I 

CFORCE 

8 

DP 

0 

DIFF3 

1 

DP 

I 

GC 

1 

DP 

I 

HC 

1 

DP 

I 

IRO 

1 

I 

I 


PS 

1 

DP 

I 

R2 

1 

DP 

I 

R3 

1 

DP 

I 


Description 

3 

3 2 

Drag force defined in PS elements 

p 3 2 " ff 3 2 
G 

H 

Flag to determine if retrograde orbit 

= -1 yes 
= .1 no 

P 

P 2 

p 3 
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ORTRAN 

ariable 

Dimension 

Type 

Input/ 

output 

Description 

UM3 

1 

DP 

I 

( a 3 + 4 ) / 2 

2 

1 

DP 

I 

a 2 

3 

1 

DP 

I 

a 3 


12 

DP 

I 

Table of averaged Fourier 
' series in 

WOL 

1 

DP 

I 

2L 

E 

1 

DP 

I 

Rotational rate of the Earth 

L 

1 

DP 

I 

P 4 

MJ 

1 

DP 

I 

V 

MUI 

1 

DP 

I 

1/U 
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C BEGIN ) 

ziz 

COMPUTE: 
CANONICAL 
FORCES 


RETURN 


Page 1 of 1 


Figure 10.- CANFOR flow chart. 
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3.3.5 CDTOJD (Subroutine) 

Purpose Compute a Julian date corresponding to a given calendar date 
Calling sequence : CALL CDTOJD ( CDATE , XJDATE ) 

Called by: INPUT 

Subroutines/functions used : None 

Named COMMON : None 

Program data : Size = 1768 C 126 iq ) words compiled 


FORTRAN 



Input/ 

variable 

Dimension 

hss. 

• output 

CDATE 

6 

DP 

I 

XJDATE 

1 

DP 

0 


Description 


A calendar date given in the form 

CDATE { 1 ) = month 
CDATE(2) = day 
CDATE(3) = year 
CDATE(4) = hours 
CDATE(5) = minutes 
CDATE(6) = seconds 

Computation range is from March 1, 1900 
through February 28, 2100 

Julian date corresponding to the given 
calendar date CDATE 
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C BEGIN ^ 

zn 


COMPUTE: 
JULIAN DATE 


RETURN 


3 


Figure 11.- CDTOJD flow chart. 
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3.3.6 COEFF ( Subroutine ) 

Purpose : Compute the binomial coefficients A( ) and the Fourier coefficients 

B( ) of the powers of cosine and sine 

Calling sequence : CALL COEFF 

Called by : POTEXP 

Subroutines/functions used : None 

Named COMMON : /EXPCOF/ A(200) ,B(200),NDEX0,NDEX(l8) ,IEXPFL 

Program data : Size = 3228 (210 iq) words compiled 


fortran 

variable 

Dimension 

Type 

Input/ 

output 

Description 

A 

200 

DP 

0 

Array containing binomial coefficients 

B 

200 

DP 

0 

Array containing Fourier coefficients 
for cosine and sine raised to a power 

IEXPFL 

1 

I 

0 

Flag to determine if the A, B arrays 
have been computed 

= 0 no 
= 1 yes 

NDEX 

18 

I 

0 

Array of pointers to A and B coefficiei 

NDEXO 

1 

X 

0 

Zero Index of table NDEX 
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C BEGIN J 

_r_ 


COMPUTE: 

BINOMIAL 
AND FOURIER 
. COEFFICIENTS 
OF THE POWERS 
OF SINE AND COSINE 


~r~ 

^ RETURN ^ 


Figure 12.- COEFF flow chart. 
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3.3*7 CONST (Subroutine) 

Purpose : Initialize the mathematical and physical constants needed to execute 

the ASOP program (refs. 13 and 14). 

Calling sequence : CALL CONST 

Called by: ASOP a , INPUT 


Subroutines/functions used : PREPS, PREPT 


Named COMMON: /CBASIC/ 

/CBODY / 
/CPRINT/ 
/DRAG / 
/END / 
/EPOCH / 
/GEO / 
/GMTROT/ 
/PERTRB/ 


PI , TWOPI ,DEG , RAD , DAY, DTOKM 

XMU , XMUI , SQTMU , SQTMUI , EPS ' 

X,I(2) ,IUNITS 

CD , AREA , XMASS , CDRAG 

STOP, ISTOP 

CDATE( 6 ) , XJDATE 

RE , C J2, CS ( 1 87 ) , SS ( 1 87 ) , IGEOFL 

WE , THETAO 

IDRAG, ILONG 


Equivalence : ( DT2 , DT 1 ) , ( THETAO , GMST ) , ( XJDATE , XJED ) 

Program data : Size = 406g (262-jq) words compiled 


FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

AREA 

1 

DP 

I 

Cross sectional surface area of 
the satellite (m^) (needed only 
if IDRAG > 1 ) 

CD 

1 

DP 

I 

coefficient of drag (needed 
only if IDRAG > 1 ) 

CDRAG 

1 

DP 

0 

A 

Drag coefficient = C rf — (m 2 /kg) 

2m 

CJ2 

1 

DP 

0 

Jg coefficient of the central 
body 

DAY 

1 

DP 

0 

Value to convert days into 
seconds, minutes-, or hours 


a 


Called only by subroutine package. 
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FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

77FM50 

Description 

DEG 

1 

DP 

0 

180/tt 

DTOKM 

1. 

DP 

0 

Converts distance into kilometers 

EPS 

1 

DP 

0 

3/2 (U J 2 Rq 2 ) .. 

GMST 

1 

DP 

I 

Initial hour angle of Earth 

IDRAG 

1 

I 

I 

Flag to determine if the drag • 
calculations are to be included 

= 0 no 

= 1 yes 

ILONG 

1 

I 

I 

Flag to determine the type of 
geopotential terms to be used 

= 0 none (two-body orbit) 

= 1 J 2 short period, and first-order 

secular terms 

= 2 Compute the mean energy due to 

geopotential terms as defined 
by NMAX and MMAX 

IUNITS 

1 • 

I 

I 

Flag to determine what units are to be 
used for the calculations 

= 1 km, sec 
= 2 nm,sec 
= 3 ft, sec 
= 4 m,sec 
= 5 km, hr 
= 6 nm,hr 
= 7 E.r. ,min 

PI 

1 

DP 

0 

ir 

RAD 

1 

DP 

0 

TT/180 

RE 

1 

DP 

0 

Equatorial radius of central body (Earth) 

SQTMU 

1 

DP 

0 

YF 

SQTMUI 

1 

DP 

0 

i/Vu 

TWOPI 

1 

DP 

0 

2TT 
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FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

WE 

1 

DP 

0 

Rotational rate of Earth 

XJED 

1 

DP 

I 

Julian date 

XMASS 

1 

DP 

I 

Total mass of the satellite (kg) 

XMU 

1 

DP 

0 

Gravitational constant for the 
central body (p) 

XMUI 

1 

DP 

0 

i/u 


48 




, Page 1 of 1 

Figure 13.- CONST flow chart. 
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1.3.8 DENSTY (Subroutine) 

urpose? Compute the atmospheric density at a given altitude above an oblate 
Earth (ref. 15) 

falling sequence : CALL DENSTY (ALT,RHO) 

■ailed by : PREPD 

ubroutines/f unctions used ? None 

amed COMMON : /CARTC / X,Y,Z,V(3) , TIME, ENERGY ,R,RI 

/CBASIC/ PI , TWOPI , DEG , RAD , DAY , DTOKM 

/D ATMOS/ FBAR , XKP , SLDAY , SADAY 

/DATMO 1 / SRAB , CRAB , SDEC , CDEC , RB , TC , TG 

./DCOEFF / A(3,3,3),B(3,9),C(3,9),D(3,4) 


rogram 

data: Size = 

672g ( 442 iq)- words compiled 

ORTRAN 



Input/ 


ariable 

Dimension 

Type 

output 

Description 


(3,3,3) 

DP 

I 

Parameters for determining the 
base altitude 

LT 

1 

DP 

I 

Altitude above an oblate Earth 


(3,9) 

DP 

I 

Parameters for determining the T OT 
density profile 


(3,9) 

DP 

I 

Parameters for computing annual 
variation 

DEC 

1 

DP 

I 

Cosine of bulge declination 

RAB 

1 

DP 

I 

Cosine of bulge right ascension 


(3,4) 

DP 

I 

Parameters for computing the seasonal 
latitudinal variation 

rOKM 

1 

DP 

I 

Converts distance into kilometers 


1 

DP 

I 

Magnitude of the position vector 
of the satellite 

3 

1 

DP 

I 

Magnitude of the diurnal change in 
the exospheric temperature 

10 

1 

DP 

0 

The computed atmospheric density 
(kg/m 3 ) 
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FORTRAN 



Input/ 


variable 

Dimension 

Type 

output 

Description 

RI 

1 

DP 

I 

Inverse R 

SADAY 

1 

DP 

I 

Magnitude of the semiannual 
density variation 

SDEC 

1 

DP 

I 

Sine of the bulge' declination 

SLDAY 

1 

DP 

I 

Magnitude of the seasonal 
latitudinal density variation 

SRAB 

1 

DP 

I 

Sine of the bulge right ascension 

TC 

1 

DP 

I 

Nighttime minimum of the global 
exospheric temperature (°K) 

TG 

1 

DP 

I 

Variation in the exospheric 





temperature due to geomagnetic 
activity (°K) 

X 

1 

DP 


Cartesian coordinates for the 

Y 

1 

DP 

1 

position of the satellite 

Z 

1 

DP 

I ) 
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BEGIN 


'kLTIl\}V£ >s ^ s N0 
> 90km 




RETURN 


Figure 14.- DENSTY flow chart. 
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3.3.9 DETERM (Subroutine) 

Purpose; Compute the first-order, long-period generating functions and its 
derivatives (ref. 8) 

Calling sequence: CALL DETERM 


Called by ; LONGPP 

Subroutines/ functions used; IL0G10 


Named COMMON: 


/DETE / SHAT ,SHATP , SHATE2 , SHATB , SHATXI ,SHATPI 
/ECC / ES.ESSQ 

/EXPCOF/ A(200) ,B(200) ,NDEXO ,NDEX(l8) 

/GEO / RE,CS2,CS(187),SS(187) 

/PSANS2/ SUM2 , SUM3 , DIFF2 , DIFF3 , GC , HC , PSSQRT , PS , QS 
/RETRO / IRO. 

/RPOOL / XO,X(l8),YO,Y(l8),PESSQO,PESSQ(9),KMMB2,KMMB21 
/S1STAV/ GIN , HOG , GPH ,BS , FS , GINSQ 
/TESS / NMAX,MMAX 
/XIPSI / XI1,PSI1 


Program data; Size = 14678 (823io) words compiled 
FORTRAN Input/ 


variable 

Dimension 

T/pe 

output Description 

A 

200 

DP 

I 

Array containing the Binomial coefficients 

B 

200 

DP 

I 

Array containing the Fourier coefficients 
for cosine and sine raised to a power 

BS 

1 

DP 

I 

b = 1 - H/G 

CS 

187 

DP 

I 

C coefficients of the geopotential model 
in the unnormalized form 

EB 

1 

DP 

0 

e VF 

ES 

1 

DP 

I 

e 

ESSQ 

1 

DP 

I 

e 2 

IESSQ 

1 

I 

*) 

Number of terms to be generated for 

IMMAX 

1 

I 

*| 

> the expansion of the Earth's geopotential 

model 

INMAX ‘ 

1 

I 

I J 

1 
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?ORTRAN Input/ 


variable 

Dimension 

Type 

output 

Description 

[RO 

1 

, I 

I 

Flag to determine if the orbit is 
retrograde 

= -1 yes 
= 1 no 

JDEX 

18 

I 

I 

Array of pointers to the A and B 
coefficients 

1MAX 

1 

I 

I 

Maximum number of zonal terms to be 
included by the geopotential model; 
maximum value of IMMAX and INMAX 

IMAX2 

1 

I 

0 

Maximum value of IESSQ; (NMAX + 1)/2 

>S 

1 

DP 

I 

P 

>SI1 

1 

DP 

I 

e sin I cos g 

IE 

1 

DP 

I 

Central body equatorial radius (R e ) 

IEOP 

1 

DP 

0 

R e /P 

5HAT 

1 

DP 

I/O 

First-order, long-period generating 
function 

5HATB 

1 

DP 

I/O > 

i 

IHATE2 

1 

DP 

I/O ] 

Derivatives of first-order, long-period 

;hatp 

1 

DP 

I/O 

/ generating function (SHAT) 

■HATPI 

1 

DP 

I/O 


iHATXI 

1 

DP 

I/O J 


:n 

1 

DP 

I 

e sin I sin g 
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C BEGIN J 

in. 

COMPUTE FIRST ORDER 
LONG PERIOD GENERATING 
FUNCTION S AND ITS 
DERIVATIVES 


A 

3S 


3S 


A 

3S 


A 

3S 


A 

3S 


3p* 3b* 3e^’ 3Y-* dx. 


1 


s i 


t 

( RETURN ^ 
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Figure 15.- DETERM flow chart. 
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3.3.10 DRAG (Subroutine) 

Purpose; Adjust the PS elements to account for the atmospheric drag pertur 
bations (refs. 6 and 7) 

Calling sequence : CALL DRAG 

Called by : PSANS 

Subroutines/functions used : None 

Named COMMON; /DRAG1 / CF0RCE(8) ,T4, TEMPO 

/PS / SIG(4),RH0(4),TAU 

Program data : Size = 115s (77 iq) words compiled 


FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

CFORCE 

8 

DP 

I 

Drag force defined in PS elements 

RHO 

4 

DP 

I/O 

PS elements Pl>P2>P3>P4 

SIG 

4 

DP 

I/O 

PS elements cr -j >P2 ,a 3 ,ff 4 

TAU 

1 

DP 

I 

Independent variable of the 
PS elements (t) 

TEMPO 

1 

DP 

I 

Second-order correction for 
density due to drag 

TLINER 

1 

DP 

I 

Change in time due to drag 

T4 

1 

DP 

I 

Magnitude of the quadratic variation 
in the mean anomaly 
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Figure 16.- DRAG flow chart. 


57 




77FM50 


3.3.11 FPRIMS. (Subroutine) 

Purpose : Compute the second-order zonal Hamiltonian and its derivatives (ref. 8) 

Calling sequence: CALL FPRIME 


Called by : LONGPP 

Subroutines /functions used : IL0G1 0 


Mamed COMMON : /ECC / ES,ESSQ 

/EXP COF/ A(200),B(200),NDEXO,NDEX(18) 

/FP / FHAT , FHATP , FHATE2 , FHATB 

/GEO / RE,CJ2,CS(187),SS(187) 

/PSANS2/ SUM2 ,SUM3 , DIFF2 , DIFF3 , GC , HC ,PSSQRT , PS , QS 

/RPOOL / YY(30) ,PESSQ0 ,PESSQ(9) j KMMB2 , KMMB2 1 

/S1STAV/ GIN , HOG , GPH ,BS , FS , GINSQ 

/TESS / NMAX,MMAX 


luivalence : 


(CS( 1 ) , c J( 1 ) ) , (SS( 1 ) ,S J( 1 ) ) 


Program data : Size = 11778 (839 -]q) words compiled 


FORTRAN 

variable 

Dimension 


Input/ 

output 

A 

200 

DP 

I 

B 

200 

DP 

I 

BS 

1 

DP 

I 

CJ 

187 

DP 

I 

ESSQ 

1 

DP 

I 

FHAT 

1 

DP 

I/O 

FHATB 

1 

DP 

I/O \ 

FHATE2 

1 

DP 

I/O > 

FHATP 

1 

DP 

I/O ) 



Array containing the binomial 
coefficients 


Array containing the Fourier 
coefficients for cosine and 
sine raised to a power 

b = 1 - H/G 


C coefficients of the 
geopotential model in the 
unnormalized . form 


Second-order zonal Hamiltonian 


Derivatives of second-order 
zonal Hamiltonian 
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FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

IESSQ 

1 

I 

r ) 

Number -of terms to be gen- 
erated for the expansion of 

INMAX 

1 

I 

I J 

the geopotential model 

NDEX 

1 

I 

I 

Array of pointers to the A 
and B coefficients ' 

NMAX 

1 

I 

I ‘ 

Maximum number of zonal terms 
to be included by the 
geopotential model; maximum 
value of INMAX 

NMAX2 

1 

I 

0 

Maximum value of IESSQ; 
(NMAX + 1 )/2 

PS 

1 

DP 

I 

P 

RE 

1 

DP 

I 

Central body equatorial ra- 
dius (R e ) 

REOP 

1 

DP 

0 

R e /p 
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Figure 17.- FPRIME flow chart. 
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3.3.12 GEOPOT (Subroutine) 

Purpose : Compute the Earth’s geopotential in a recursive manner (ref. 16) 

Calling sequence : CALL GEOPOT (POT) 

Called by : PSTOX, XTOPS 


Subroutines/functions used: None 


Named COMMON: /CARTC / 

X,Y,Z,VX(3),TIME,ENERGY,R,RI 

/CBASIC/ ' 

PI , TWOPI , DEG , RAD ,DAY , DTOKM 

/ CBODY / 

XMO , XMUI , SQTMU , SQTMUI , EPS 

/GEO / 

RE , C J2 , CS( 1 87 ) , SS ( 1 87 ) , IGEOFL 

/GMTROT/ 

WE , THETAO 

/PERTRB/ 

IDRAG,ILONG 

/RPOOL / 

POO,PO(19) ,P10 ,P1 (19) ,P20,P2(19) ,CTIL0,CTIL(19) , 
STILO , STIL ( 1 9 ) , R2I , R3I 

/TESS / 

NMAX .MMAX ' 


Equivalence : (CTIL( 1 ) ,CTIL1,) , (STIL( 1 ) ,STIL1 ) 


Program data : Size = ^30g (280jg) words compiled 


FORTRAN 

variable 

Dimension 

M 

Input/ 

output 

Description 

CS 

187 

DP 

I 

C coefficient of the 
geopotential model in the 
unnormalized form 

EPS 

1 

DP 

I 

= 3/2 (V J 2 R e 2 ) 

ILONG 

1 

I 

I 

Flag to determine the type 


of .geopotential terms to be 
included 


= 0 none (two-body orbit) 

= 1 J ; 2 short period, and 
first-order secular 
terms 

= 2 Compute the mean energy 
due to the geopotential 
terms as defined by NMAX 
and MMAX 
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FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

MMAX 

1 

I 

I 

Total number of tesseral 
terms, to be included -by -the 
geopotential model (value is 
needed only when ILONG =2) 

POT 

1 

DP 

0 

Magnitude of the Earth's 
gravitational potential 

R 

1 

DP 

I 

Magnitude of the position 
vector of the satellite 

RE 

1 

DP 

I 

Equatorial radius - of the 
central body (Earth) 

RI 

1 

DP 

I 

Inverse of R 

ss ■ 

187 

DP 

•I 

S coefficients of the 
geopotential model in 
the unnormalized form 

THETAO 

1 

DP 

I 

Initial hour angle of the 
Earth 

TIME 

1 

DP 

I 

Elapsed time 

TWOPI 

1 

DP 

I 

2ir 

WE 

1 

DP 

I 

Rotational rate of the Earth 

X 

1 

DP 

I 

X-Component of the Earth 
inertial position vector 

XMU 

1 

DP 

I 

U 

Y 

1 

DP 

I 

Y-Component of the Earth 
inertial position vector 

Z 

1 

DP 

I 

Z -Component of the Earth 
■inertial position vector 
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3.3.13 ILOGIO (Subroutine) 

Purpose' : Determine the number of terms to be included in the expansion for the 

Earth's geopotential model 

Galling sequence : CALL ILOGIO (X, MAX, IMAX) 

Called by : DETERM, FPRIME,POTEXP 

Subroutines/functions used: none 


Named COMMON : none 

Program data : Size = 106g (70 -jq) words compiled 


FORTRAN 

variable Dimension Type 
IMAX 1 I 

MAX 1 I 

X 1 DP 


Input/ 

output . Description 

0 Number of terms , between 1 

and MAX, to be included in 
the expansion for Earth's 
geopotential model 

1 Maximum value of IMAX 

I Small parameter used in 

determining number of terms 
to be’ included in the 
expansion 
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3.3.1# INITAL 

“ 7 T ~ ' U" 

Purpose : In^ti^lizp the coefficients for the Jacchia 71/Lineberry atmospheric 

^ r *' density model 


Calling sequence : CALL INITAL (TIME) 

Called by : flR$PD 


Subroutines/functions used : SACT , SUN 

"* T,.- y / - - -] ' t - 'I - -t 1 


Named COMMON ; /CBASIC/ 
' ■’ /DATfTOS/ 

/DATM01/ 
/DCOEFF/ 
/EPQCH / 
/SUNPOS/ 


PI , TWOPI ,PEG , RAP,DAYSEC f DTOKM 

FBAR . AKP , SLDAY , SADAY 

SRAB , CRAB , SDEC , CDEC , RB , TC , TG 

A(3,?,3),B(3,9),C(3,9) ? D(3,#) 

CDATEC6) iXJPATE 

)f1S,X?^»X3S,RS,RAS,DECS 


Program data: 

T T T r» '' . 


Size = 255$ (173 iq) word? compiled 


FORTRAN 

variable 

Dimension 

type 

Input/ 

output 

Description 

1 P 

A 

(3,3,3) 

DP 

0 

Parameters for determining 
the base altitude 

i\KP 

1 

DP 

I 

Averaged value for the 
geomagnetic index 

B 

(3,9) 

DP 

0 

Parameters for determining 
the To, density profile 

9 

(3,9) 

DP 

0 

Parameters for computing an- 
nual variation 

CDEC 

1 

DP 

0 

Cosine of bulge declination 

CRAB 

1 

DP 

0 

Cosine of bulge right 
ascension 


(3,#) 

DP 

0 

Parameters for computing 
the seasonal latitudinal 
variation . 

DAY 

1 

DP 

0 

Julian day number for 
which the solar activity 
coefficients are desired 

DAYSEC 

1 

DP 

I 

Converts days into seconds, 
minutes, or hours 
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FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

FBAR 

1 

DP 

I ‘ 

Averaged value for__the solar 
flux coefficient, F-|q .7 

RAD 

1 

DP . 

I V 

180 n 

HAS 

1 

DP 

I 

| Right ascension of the Sun 

RB 

1 

DP 

• . .0 

] Magnitude of the diurnal 
change in the exospheric 
temperature 

RS 

1 

DP 

I 

Magnitude of position vector 
of the Sun 

SADAY 

1 

DP 

0 

Magnitude of the semiannual 
density variation 

SDEC 

1 

DP 

0 

Sine of bulge declination 

SLDAY 

1 

DP 

, 0 

Magnitude of the seasonal 
latitudinal density variation 

SRAB 

1 

DP 

0 '• 

t 

Sine of bulge right ascension 

TC 

1 

DP 

0 

Night time minimum of the 
global exospheric temperature 
(°K) 

TG 

1 

DP 

0 

Variation in the exospheric 
temperature due to geomagnetic 
activity (°K)- 

TIME 

1 

DP 

I 

Elapsed time of epoch 

TWOPI ' 

1 

DP 

I 

2TT 

XJDATE 

1 

DP 

I 

Julian date of epoch 

X3S 

1 

DP 

I 

Position vetor of the Sun 
in the Earth's inertial 
equatorial system (X1S, 
X2S, X3S) 
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Figure 20,- INITAL flow chart. 
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3.3.15 INPUT (Subroutine) 

Purpos e: Read the input data from the NAMELIST statement, set the default 

values, and initialize all required COMMON block variables. 


Calling sequence: CALL INPUT ($20.) 

4 


MAIN 


Subroutines/functions used: AEIXYZ, CDTOJD, CONST, OUTPUT, XTOPS 


Named COMMON: /CARTC 

/ 

X( 6), TIME, ENERGY, R,RI 

/CBASIC 

/ 

PI , TWOPI ,DEG, RAD ,DAY ,DT0KM 

/CPRINT 

/ 

PRINT , IPRINT , IPSPRT , IUNITS 

/DRAG 

/ 

CD, AREA, XMASS, DRDR 

/END 

/ 

STOP, ISTOP 

/EPOCH 

/ 

DATE(6),XJDATE 

/KEPLER 

/ 

EL(6) 

/PERTURB/ 

IDRAG, ILONG 

/TESS 

/ 

NMAX, MM AX 


NAMELIST statements: /INPUT/ EL, IEL, STOP, ISTOP, PRINT, IPRINT, DATE, IDRAG, CD, 

AREA , XMASS , ILONG, NMAX ,MMAX , IPSPRT , IUNITS 

Program data : Size = 4070 ( 263 -]q) words compiled 


FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

AREA 

1 

DP 

I 

Cross sectional surface area 
of the satellite (needed only 
if IDRAG > 1) (m 2 ) 

CD 

1 

DP 

I 

Cbefficient of drag (needed 
only if IDRAG > 1 ) 

DATE 

6 

DP 

I 

Calendar date of epoch 


DATE( 1 ) 

= month 

(2) 

= day number 

(3) 

= year 

(4) 

= hour 

(5) 

= minutes 

(6) 

= seconds 
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FORTRAN 



Input/ 

variable 

Dimension 

T ype 

output 

EL 

* 

DP 

I/O 


ENERGY ) DP . 0 

?DRAq 1 X I 


IEJL 


I 


I 


3ILQNq 1 I I 


Description 

Initial conditions of the sat- 
ellite given in Keplerian ele- 
ments of Cartesian coordinates 
on output it will contain the 
Keplerian elements 

EL(1) a X or a 

(2) = Y or e 

(3) = Z or i 

(4) a x or w 

(5) = Y or ft 

(6) = Z or M 

Total energy of the satel- 
lite; initially set to 0 

Flag to determine if the drag 
calculations are to be 
included 

= 0 no 
a 1 yes 

Flag to determine if the 
input values of EL are given 
as Keplerian elements or 
Cartesian coordinates 

= 1 Keplerian 
a 2 Cartesian 

Flag to determine type of. 
geopotential terms to ' 
be included 

= 0 none (two -body) 

= 1 J 2 short period, and 

first-order secular terms 
= 2 Compute the mean energy 
due to the geopotential 
terms as defined by NMAX 
and MAX 
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FORTRAN 

variabl e Dimension Type 
IPRINT 1 I 


IPSPRT 1 I 


ISTOP 1 I 


IUNITS 1 I 


MM AX 1 I 


NMAX 1 I 

PRINT 1 DP 


Input/ 

output - Description - 

I Flag to determine if the in- 

termediate printout is to be 
done at a PRINT value of days 
or revolutions 

= 0 no intermediate printout . , 
= 1 days 
= 2 revolutions 

I Flag to determine if the PS 

elements are to be included 
with all output 

= 0 no 
= 1 yes 

I Flag to determine if the STOP 

condition is days or revolu- 
tions 

= 1 days 
= 2 revolutions 

I Flag to determine what calcu- 

lation constants are to be 
used 

- 1 km, sec 
= 2 nm,see 
= 3 ft, sec 
= 4 m,sec 
= 5 km, hr 
= 6 nm,hr- 
.= 7 E.r.,min 

I Maximum number of tesseral 

terms to be included (needed 
only if ILONG > 2) 

I Maximum number of zonal terms 

to be included (needed only 
if ILONG > 2) 

I • Increment at which the inter- 

mediate printout is desired 
(needed only if IPRINT > 1) 
(days or revs') 
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FORTRAN 



Input/ 


variable 

—,— 1 - 

Dimension 

Type 

output 

Description 

RAD 

1 

DP 

I 

tr /180 

£TQP 

1 

DP 

I 

Final stop value at which 

’v 



output is desired (days or 
revs) 


TIME 

1 

DP 

0 

Physical time; initially set 
to 0 (hrs or min or sec) 

X 

6 

DP 

0 

Initial Cartesian state 
vector 


- 



X(1) + X(3) = X 





X(A) -> X(6) = V 

xjdate 

1 

DP 

I 

Julian date 

XMAS5 

1 

DP 

I 

Initial mass- of the satellite 
(kg) 
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Q BEGIN ^ 


SET DEFAULTS 


READ $INPUT 
PRINT $ INPUT | 
CALL CDTOJD 
CALL CONST 



DEG. TO RAD. 
CALL AEIXYZ 

1 

□ 

* 

f 


TIME = 0 
ENERGY ■ 0 



CALL XTOPS 
CALL OUTPUT 


a 

SOLICIT FOR 
CHANGES TO 
INPUT 


INPUT 

*X' 5 *C', or 'S' 



Page 1 


Figure 21.- INPUT flow chart. 
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'3. 3. If JPTOCD (Subroutine) 

, ‘ * - | *r r ■ " M ' 1 ■ — 

Purpose ; Determine the calendar date corresponding to a given Julian date 
C&lling sequence : CALL JDTOCD (CD ATE, X JD ATE ) 

Ca lled by : OUTPUT 

1 ' j ' 

Spb routines/functions used: None 

-■h t — f 7 

Named COMMON : None 

Program data : Size = I60g (112 10 ) words compiled 


FORTRAN 

variable Dimension 

pDATE 6 


XJDATp 1 


Type 

DP 


Input/ 

output 


DP 


Description 


The calendar date 
corresponding to XJDATE, 
given in the form 


CDATE ( 1 ) 

= month 

(2) 

= day 

( 3 ) 

= year 

( 4 ) 

= hours 

(5) 

= minutes 

( 6 ) 

= seconds 

A Julian 

date 
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C BEGIN J 


COMPUTE: 

CALENDAR 

DATE 


RETURN 


Figure 22.- JDTOCD flow chart. 


Page 1 of J 
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3-3.17 LQNGPP (Subroutine) 


Purpose ': Cqmppte the first-order, zonal, long periodic perturbations and 

second-order zonal, secular perturbations (ref. 9) 


Calling, sequence : 
Called by ; PSANS 


CALL LONGPP-(NN) 


Subroutines/f unctions used: DETERM, FPRIME 


Named COMMON: 


/CBODY / 

XMU , XMUI , SQTMU , SQTMUI , EPS - 

/CONSTtf/ 

TW03 ,BY3 >BY6 , CN 

/DENS / 

B( 10) ,DS1 ,DS2,DC1 ,DC2,00,OC2,OS2,OC1 ,OS1 

/DETE / 

SHAT ,SHATP ,SHATE2 ,SHATB ,SHATXI ,SHATPI 

/ECC / 

ES,ESSQ 

/FP / 

FHAT,FHATP ,FHATE2 ,FHATB 

/HAMDS / 

DSF(4) ,DSB(4) 

/PS / 

SIG(4),RHO(4),TAU 

/PSANSV / 

FACTOR (4), SIGINK8) , 

/PSANS 1/ 

ETA1 , ZETA1 , TWOL , IQL , FAK 

/PSANS2/ 

SUM2 , SUM3 , DIFF2 , DIFF3 ,G,H, PSSQRT ,PS , QS 

/PSANS 3/ 

QC 

/RP00L / 

XX(30) 

/S1STAD/ 

GC(8) ,P(8) ,Q(8) ,HC(8) ,QCV(8) 

/S1STAV/ 

GIN , HOG , GPH ,BS , FS , GINSQ 

/TESS / 

NMAX,MMAX 

/XIPSI / 

XII ,PSI1 


Equivalence : (L,RH0(4) ) , (DSB2( 1 ) ,XX( 1 ) ) , (DSD( 1 ) ,XX(5) ) , (DSP( 1 ) ,XX(9) ) , (DSE( 1 ) , 

XX( 13) ) , (DSQ( 1 ) ,XX( 17) ) , (DSDEL( 1 ) ,XX(21 ) ) ,.(DSDEL1 ( 1 ) ,XX(25) ) 

Program data : Size = 17058 (965-jg) words compiled 

A description of the mathematical symbols used and their' relationship to one an- 
other i? given in references 8 and 17; a brief description can also be found 
in Appendices E and F. Therefore, only a brief mathematical description will be 


given. 

FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

?s 

1 

DP 

I 

b = 1 - H/G 

CN 

1 

DP 

I 

+1 depending on value of NN 

DSB 

4 

DP 

I 

9b/38 k k = 1, 2, 3, 4 

DSF 

4 

DP 

I 

3f/98 k k = 1, 2, 3, 4 


i 
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FORTRAN Input/ 


variable 

Dimension 

Type 

output 

Description 

EPS 

1 

DP 

I 

E = 3/2 (VJ 2 R e 2 ) 

ES 

1 

DP 

0 

e 

ESSQ 

1 

DP 

0 

e 2 . 

FACTOR 

4 

DP 

0 

Derivatives of the DS 
Hamiltonian and its combina- 
tions (A<|, A 2 , Aj and Ajj in 
Appendix F) 

FAK 

1 

DP 

I 

<2L)-3 /2 

FHAT 

1 

DP 

I/O 

Second-order zonal 
Hamiltonian 

FHATB 

1 

DP 

I/O > 

1 Derivatives of second-order 

FHATE2 

1 

DP 

I/O ' 

^ zonal Hamiltonian 

FHATP 

1 

DP 

I/O J 

• 

FS 

1 

DP 

I 

f = 1/pq 

GC 

8 

DP 

I 

9G/aCT k , 9G/3p k k = 1, 2, 3, 4 

GIN 

1 

DP 

I 

G- 1 

GINSQ 

1 

DP 

I 

G~ 2 

GPH 

1 

DP 

I 

G + H 

HC 

8 

DP 

I 

9H/3a k , 3H/3p k 
k = 1, 2, 3, 4 

HOG 

1 

DP 

I 

H/G 

L 

1 

DP 

I/O 

L = Pl| = CTg 

NMAX 

1 

I 

I 

Maximum number of zonal terms 
to be included 

NN 

1 

I 

I 

, Flag determining if 

initializing or computing . 

= 0 initializing 
= 1 computing 
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FORTRAN 

variable 

Dimension 

'type 

Input/ 

output 

Description 

OC1 

1 

DP 

M 

Coefficients of Fourier 
series describing density 

0S1 

1 

DP 

0 ) 

due to long-period J2 
changes in height 3 

P 

8 

DP 

I 

3p/3(5 k , 3p/3p k 
k = 1, 2, 3, 4 

PS 

1 

DP 

I 

P 

PSI1 

1 

DP 

0 

t!>1 = e sin I cos g 

PSSQRT 

1 

. DP 

I 

W 

Q 

8 

DP 

I 

3q/3a k> 3q/3p k 
k = 1, 2, 3, 4 

QC 

1 

DP 

I 

Q 

QCV 

8 

DP 

I 

3Q/3a k , 3Q/3p k 
k = 1, 2, 3, 4 

QS 

1 

DP 

I 

Q 

RHO. 

4 

DP 

I/O 

P-j , . . . , P4 (see L) 

SHAT 

1 

DP 

I/O 

First-order , long-period 
generating function 

SHATB 

1 

DP 

I/O ' 

V 

SHATE2 

1 

DP 

I/O 

j Derivatives of first-order 

SHATP 

1 

DP 

I/O 

> long-period generating 

function 

SHATPI 

1 

DP 

I/O 


SHATXI 

1 

DP 

I/O , 


SIG 

4 

DP 

I/O 

<*1f " * ) CT]| 

SQTMU 

1 

DP 

I 

Vu 

SQTMUI 

1 

DP 

I 

i/ViT 

TW03 

1 

DP 

I 

2/3 


^Reference 7, pp. 18-22. 
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FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

TWOL 

1 

DP 

I 

2L 

XII 

1 

DP 

0 

X-j = e sin I sin g 

XMU 

1 

DP 

I 

V . 

XMUI 

1 

DP 

I 

1/y 
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COMPUTE LONS PERIOD 


GENERATION FUNCTION S } * 

as * as * 

AND PS DERIVATIVES — and — 1 
3cr 3p 




COMPUTE DS 
DERIVATIVES 



COMPUTE SECOND ORDER 
HAMILTONIAN AND 
ITS DERIVATIVES 


RETURN 


ADD LONG PERIOD 
CONTROBUTION TO 
a AND p 


( RETURN ^ 


Figure 23.- LONGPP flow chart. 


80 







77FM50 


3.3.18 MATIN (Subroutine) 

Purpose : Invert an n x n matrix and/or solve the matrix equation Ax = B 

Calling sequence : CALL MATIN (A, N,B,M,KEY, DETERM) 

Called- by : PREPD 

Subroutines/functions used : - None 
Named COMMON : None 

Program data : Size = 573s (379-|q) words compiled 


FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

A 

(10,10) 

DP 

I 

Matrix to be inverted 

B 

10 

DP 

I 

Column matrix to be 
multiplied by A 

DETERM 

1 

DP 

0 

Determinant 

KEY 

1 

I • 

0 

Flag for singular matrix 

M 

1 

I 

I 

Flag to determine if the B 
matrix is to be used 

= 0 no 
= 1 yes 

N 

1 

I 

I 

Size of A matrix and length 
of B 


81 



BEGIN 


INVERT MATRIX 
COMPUTE DETERMINANT 
SOLVE Ax = B 


C 


RETURN 


Figure 24.- MATIN flow chart. 


Page 1 of 1 
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3.3.19 MTOECC (Subroutine) 

Purpose : Convert mean anomaly to eccentric anomaly and compute its sine and 

cosine 

Calling sequence : CALL MTOECC (XM,E,EA,SINEA,COSEA) 

Called by : AEIXYZ, SDN 

Subroutines/functions used : None 

Named COMMON : None 

Program data : Size = I63g (115 -iq) words compiled 


FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description ' 

COSEA 

1 

DP 

0 

Cosine of the eccentric 
anomaly (cos E) 

E 

1 

DP 

I 

Orbital eccentricity (e) 

EA 

1 

DP 

0 

Eccentric anomaly (rad) 

SINEA 

1 

DP 

0 

Sine of the eccentric anomaly 
(sin E) 

XM 

1 

DP 

I 

Mean anomaly (M)(rad) 
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Figure 25.- MTOECC flow chart. 
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3.3.20 OUTPUT (Subroutine) 

Purpos e: Print all desired output during the execution of the ASOP program; it 

contains all output formats and does all unit conversions required for 
output . 

Calling sequence : CALL OUTPUT (IFORM) 


Called by : INPUT, MAIN 

Sub routines /functions used : JDTOCD, XYZAEI 


Named COMMON: /CARTC / 

/ CBASIC/ 
/CP HINT/ 
/END / 
/EPOCH / 
/KEPLER/ 
/PS / 


X(6) , TIME, ENERGY, R,RI 

PI , T WOPI , DEG , RAD , D AY , DTOKM 

PRINT , IP , IELPRT , IU 

STOP, ISTOP 

CDATE( 6 ) , XJD ATE 

XKEP(6) 

SIG( 4 ) , RHO ( U ) , TAU , TAUMAX , TAUINT 


Equivalence : (SEC,CDATE(6)) 


Program data : Size = 7108 (^56io) words compiled. 


FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

CDATE 

6 

DP 

I 

Calendar date of output in 
the form 

CDATE(I) = day number 

(2) = month 

(3) = year 

(4) = hours 

(5) = minutes 

(6) = seconds 

DAY 

1 

DP 

I 

Value to convert days into 
hours, minutes, or seconds 

DEG 

1 

DP 

I 

180/tt 

ENERGY 

1 

DP 

I 

Total energy of the physical 
system 

IELPRT 

1 

I 

I 

Flag to determine if the PS 
elements are to be printed 

= 0 no 
= 1 yes 


85 



77 FM 50 


FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

IFORM 

1 

I 

I 

1 

Flag to determine if the ini' 
tial or final condition mes- 
sage is to be printed. 

= 1 initial condition raessag< 
= 2 no message (intermediate 
print ) 

= 3 final condition message 
or reentry condition 
message 

IP 

1 

I 

1 

Flag to determine what print 
condition is being used 

= 1 days 
= 2 revolutions 
= 3 reentry 

ISTOP 

1 

I 

I 

Flag to determine the final 
stop condition 

= 1 days 
= 2 revolutions 
= 3 reentry 

10 

1 

I 

I 

Pointer to the VEL and DST 
character arrays 

.RHO 

4 

DP 

I 

PS elements p-j -* P4 

RH 0 ( 1 ) = p-, RH 0 ( 3 ) = P3 

RH 0 ( 2 ) = p 2 RH 0 ( 4 ) = P4 

SIG 

4 

DP 

I 

PS elements a-j -*• (J4 

SIG( 1 ) = a- 1 SIG( 3 ) = a 3 

SIG( 2 ) = a 2 SIG( 4 ) = 04 

TAU 

1 

DP 

I 

■Independent variable of the 
PS elements set (rad) 

TAUINT 

1 

DP 

I/O 

Initial value of TAU for 
which the force model is 
valid. Initially set to 0 . ' 

t;ime 

1 

DP 

I 

Physical time 

TWOPI 

1 

DP 

I 

2it 
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FORTRAN 

variable 

Dimension 

J 

j’ype 

Input/ 

output 

77FM50 

Description 

X 

6 

DP 

I 

Cartesian state vector 
X(1) + X(3) = X 
X(4) ->■ X(6) = V 

XJDATE 

1 

DP 

I 

Julian date 

XKEP 

6 

DP 

I/O 

Keplerian elements 


XKEP(1) = a 

(2) = e 

(3) = i 

(4) = b) 

(5) = ft 

(6) = M 
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(RETURN) 

Figure 26.- OUTPUT flow chart. 
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3.3.21 POTS XT (Subroutine) 

Purjjose: Compute the mean energy due to tesseral and sectorial geopotential har- 

monics (ref. 10 ) 

Calling sequence: CALL POTEXP 


Called by; ASOP 


Subroutines/functions used : COEFF , IL0G1 0 , RECUR , TIMEXP 


Named COMMON : /CBASIC/ 

/CBODY / 
/ EXPCOF/ 
/GEO / 
/GMTROT/ 
/PS / 
/PSANS 1 / 
/PSANS2/ 
/PSANS3/ 
/RETRO / 
/RPOOLA/ 


/TESS / 


PI , TWOPI , DEG , RAD , DAY , DTOKM 

XMU , XMUI , SQTMU , SQTMUI , EPS 

A( 200 ) , B ( 200 ) , NDEXO , NDEX (18), IEXPFL ■ 

RE , C J2 , CC ( 1 87 ) ,SS( 1 87 ) , IGEOFL 

WE , THETAO 

SIG(4) ,RHO( 4) ,TAU 

SSIG1 , CSIG1 , TWOL ,XXQL ,FAK 

SUM2 ,SUM3 , DIFF2 ,DIFF3 , GC ,HC ,PSSQRT ,PS , QS 

QCAP 

IRO 

FNO(24),FN(19,19),FN10(20),FN1(19,19),FM20(20), 
FN2 ( 1 9 , 1 9 ) , BET AO ,BETA( 18), DELTAO , DELTA( 18), 
ETA(20) ,ZETA0 ,ZETA(20) ,ALPH0,ALPH(36) , GAMMO , GA 
XI0,XI(38) ,PSX0,PSI(38) ,GCAP0,GCAP(17) 

NMAX,MMAX 


Program data : Size = 22140 (1164-jq) words compiled 


FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

A 

200 

DP 

I 

Array containing binomial 
coefficients 

ALPH 

37 

DP 

I/O 

Array generated recursively, 
where ALPH(1) = P P 3 

B 

200 

DP 

I 

Array containing Fourier 
coefficients 

BETA 

19 

DP 

I/O 

Array generated recursively, 
, where BETA(1) = 

cos (WE • + THETAO) 

BS 

1 

DP 

0 

Vi - (H/G ) 2 

CC 

187 

DP 

I 

C coefficients of the geo- 
potential model in the 


unnormalized form 
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FORTRAN 

variable 

Dimension 


Input/ 

output 

Description 

CSIG 1 

1 

DP 

I 

COS O-j 

DJSLTA 

19 

DP 

I/O 

Array generated recursively, 
where DELTA ( 1 ) = 
sin (WE • 04 + THETAO) 

ECC 

1 

DP 

0 

e 

ECC 2 

1 

DP 

0 

e2 

ETA 

21 

DP 

I/O 

Array generated recursively, 
where ETA( 1 ) = -Q ^2 

GAMM 

37 

DP 

I/O 

Array generated recursively 
where GAMM( 1 ) = -P 03 

GC 

1 

DP 

I 

G 

HC 

1 

DP 

I 

H 

IEXPFL 

1 

I 

I 

Flag to determine if A and B 
arrays have been computed 

= 0 no 
= 1 yes 

INMAX 

1 

I 

I/o' 

Number of terms to be 
generated in BETA, DELTA and 
XM 

IQMAX 

1 

I 

I/O 

Number of terms to be 
generated in ALPH and GAMM 

IRMAX 

1 

I 

I/O 

Number of terms to be 
generated in ZETA and ETA 

IRO 

1 

I 

I 

Flag to determine if the 
orbit is retrograde 

= -1 yes 
= 1 no 

ISMAX 

1 

I 

0 

Number of terms to be 
generated in PSI and XI 

MAX 

1 . 

I 

0 

Maximum value of IRMAX 
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FORTRAN Input/ 

variable Dimension Type output Description 


MMAX 

1 

I 

I 

Maximum number of tesseral 
terms to be included by the 
geopotential model; maximum 
value of INMAX; 2*MMAX is 
maximum value of IQMAX 
(Needed only when ILONG =2) 

NDEX 

18 

I 

I 

Array of pointers to the A 
and B coefficients 

NDEXO 

1 

I 

I 

Zero Index to the NDEX array 

PS 

1 

DP 

I 

P 

PSI 

39 

DP 

I/O 

Array generated recursively, 
where PSI(1) = cos cr -j 

QCAP 

1 

DP 

I 

Q 

QS 

1 

DP 

I 

q 

RE 

1 

DP 

I 

Equatorial radius of the cen- 
tral body (Re) 

REOP 

1 

DP 

0 

Rg/p 

RHO 

4 

DP 

I/O 

PS elements Pi,P 2 >P3»P4 
Note P4 = energy 

SIG 

4 

DP 

I 

PS elements o-j 

SS 

187 

DP 

I 

S coefficients of the 
geopotential model in 
the unnormalized form 

SSIG1 

1 

DP 

I 

sin cr^ 

SUM2 

1 

DP 

I 

l(a 2 2 + P2 2 ) 

THETAO 

1 

DP 

I 

Initial hour angle of the 
Earth 

WE 

1 

DP 

I 

Rotational rate of the Earth 

XI 

39 ‘ 

DP 

I/O 

Array generated recursively, 
where XI( 1 ) = sin a ^ 
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FORTRAN 



Input/ 


variable 

Dimension 

Type 

output 

Description 

XIQL 

1 

DP 

I 

U/V2IT 

XM 

(7,18) 

' DP 

0 

Fourier coefficients of the 
time expansion 

XMU 

1 

DP 

I 

U 

XNU 

1 

DP 

0 

Ratio of the frequency of ro- 
tation of the satellite to 
the rotation rate of the 
Earth 

ZETA 

21 

• DP 

I/O 

Array generated recursively, 
where ZETA( 1 ) = Q P2 
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Figure 27.- POTEXP flow chart. 
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3.3.22 PRBPD (Subroutine) 

Purpose : Initialize the ASOP atmospheric density model coefficients 

(refs. 6 and 7) 

Calling sequence : CALL PREPD 

Called by: ASOP 


Subroutines/functions used: CANFOR , DENSTY , INITAL , MATIN } PSTOX , TABLE 


Named COMMON : /CARTC / 

/CBASIC/ 
/CBODY / 
/ DBETAS/ 
/DENS / 
/DRAG / 
/DRAG1 / 
/END / 
/ FORSAV/ 
/GEO / 
/GMTROT/ 
/PERTRB/ 
/PS / 
/PSANSV/ 
/PS ANSI/ 
/PSANS2/ 
/PSANS3/ 
/RETRO / 
/SUNPAR/ 
/SUNPOS/ 


X( 6), TIME, ENERGY, R,RI 

PI , TWOPI , DEG , RAD , DAY , DTOKM 

XMU , XMUI , SQTMU , SQTMUI , EPS 

BETA1 ,BETA2,BETA3,BETA4 

B( 10) ,DS1 ,DS2,DC1 ,DC2,00,OC2,0S2,0C1 ,OS1 

CD , AREA , XMASS , CDRAG 

CFORCE( 8 ) , T4 , TTL , TEMPO 

STOP , ISTOP 

Z(6),ZC(9,6),ZS(9,6),DZ(6),DZC(9,6),DZS(9,6) 
RE,CJ2,CS( 187) ,SS( 187) , IGEOFL 
WE , THETAO 
IDRAG , ILONG 

SIG(4),RH0(4),TAU,TAUMAX 
FACT0R(4) ,SIGINI(8) 

SIN1 ,COS1 , TWOL , XIQL , FAK 

S0M2 ,SUM3 , DIFF2 ,DIFF3,G ,H ,PSSQRT ,PS , QS 

QC 

IRO 

XNS 

XS , YS ,ZS , RS , R AS , DECS 


Program data : Size = 221 2g ( 1 1 62 i ) words compiled 


FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

• Description 

ALT 

1 

DP 

0 

Calibration altitude for the 
density model 

ALTO 

1 

DP 

0 

Reference altitude for the 
density model 

B 

1 

DP 

0 

Coefficient of the density 
model (ref. 7) 

BETA1 

1 

DP 

0 

8 

BETA2 

1 

DP 

0 

8 2 

BETA3 

1 

DP 

0 

8 3 
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FORTRAN Input/ 


variable 

BETA4 

Dimension 

Type 

DP 

output 

Description 

i 

0 

B 4 


CDRAG 

i 

DP 

I 

Drag coefficient 

CFORCE 

8 

DP 

I/O 

Drag force defined in PS 
elements 

COS1 

1 

DP 

. I/O 

cos < 

? 1 

DC1 

1 

DP 

°) 

Coefficients of the Fourier 
> series describing the 

DC2 

1 

DP 

0 ] 

l diurnal bulge (ref . 7 ) 

DECS 

1 

DP 

I 

Declination of the Sun 

DIFF3 

1 

DP 

I 

P 3 2 

-«3 2 

DS1 

1 

DP 

° ; 

| Coefficient's of the Fourier 

> series describing the 

DS2 

1 

DP 

0 

) diurnal bulge (ref. 7) 

DTOKM 

1 

DP 

I 

Converts distance into 
kilometers 

DZ 

6 

DP 

I 

*} 

1 Coefficients of Fourier 

series expansion of 

DZC 

(9,6) 

DP 

I 

\\ 

> the disturbing 
i function in PS 

DZS 

(9,6) 

DP 

I 


) elements (ref. 7) 

FACTOR 

4 

DP 

I 

Derivatives of DS Hamiltonian 
and its combinations (Ai,A 2 , 
Ag and Aij in appendix F) 

FAK 

1 

DP 

I 

(2L)-3 /2 

G 

1 

DP 

I 

G 


H 

1 

DP 

I 

.H 


IDRAG 

1 

DP 

0 

Flag to determine if drag cal- 
culations are to be included 

= 0 no 
= 1 yes 
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FORTRAN Input/ 


variable 

Dimension 

Type 

output 

Description 

IRO 

1 

DP 

I 

Flag to determine if the 


orbit, is retrograde- 


= -1 yes 
= 1 no 


0C1 

1 

DP 

I/O ^ 


3C2 

1 

DP 

1/0 

Coefficients of the Fourier 
series describing density 

DS1 

1 

DP 

1/0 ! 

variation due to Jg changes 
I in height (ref. 7 ) 

DS2 

1 

. DP 

I/O 


)0 

1 

DP 

I ) 


PSSQRT 

1 

DP 

I 

P 

5C 

1 

DP 

I 

Q 

3S 

1 

DP 

I 

<3 

i 

1 

DP 

I 

Magnitude of the position 
vector of the satellite 

iAD 

1 

DP 

I 

180 /IT 

IAS / 

1 

DP 

.1 

Right ascension of the Sun 

\E 

1 

DP 

I 

Equatorial radius of the cen- 
tral body 

IHO 

4 

DP 

I 

PS elements' Pi,p 2 ,p 3 ,P 4 

moDi 

1 

DP 

0 

Density at the calibrated 
altitude 

IH0D2 

1 

DP 

0 

Density at the reference 
altitude 

>16 

4 

DP 

I/O 

PS elements ^, 02 , 03,04 

!IN1 

'1 

DP 

I/O 

sin ct-j 

►QTMU • 

1 

DP 

I 

VTT 

►TOP 

1 

DP 

I 

Final stop value at which 
output is desired 
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FORTRAN 

variable 

Dimension 

T ype 

Input/ 

output 


Description 

SUM3 

i 

DP 

I 

^(<*3^ + P3^) 

TAUMAX 

i 

DP 

0 

Range. of validity for the 
force model 

TEMPO ‘ 

i 

DP 

0 

Second order correction for 
density due to drag 

TIME 

i 

DP 

I 

Elapsed time of epoch 

TTL 


. DP 

0 

Change in time due to 
drag (At) 

TWOL 

i 

DP 

i 

2L 


TWOPI 

i 

DP 

i 

2it 

■ 

T4 

i 

DP 

0 

Magnitude of the quadratic 
variation in the mean anomaly 

WE 

i 

DP 

i 

Rotational rate of the Earth 

XIQL 

i 

DP 

i 

U/ 2L 

XMU 

i 

DP 

i 

Gravitational constant of the 
central body (]i) 

XNS 

i 

DP 

i 

Mean 

motion of Sun 

Z 

6 

DP 

i 

1] 

| ' Coefficients of Fourier 
series expansion of 

ZC 

(9,6) 

DP 

i 

XJ 1 

► the disturbing function 
in PS elements (ref. 7) 

ZS 

(9,6) 

DP 

i 
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3.3.23 PREPS (Subroutine) 

Purpose : Establish the parameters needed to calculate the position of the Sun 

(ref. 18) 

Calling sequence : CALL PREPS 

Called by : CONST 

Subroutines/functions used : none 

Named COMMON : /CBASIC/ PI , TWOPI , DEG , RAD ,DAY ,DTOKM 

/EPOCH / CDATE( 6 ) , X JDATE 

/SUNPAR/ XNS,XLS0,A,E,AER00T,B1(2),B2(2),B3(2) 

Equivalence : (C0,B1(1)) 

Program data : Size = 276q ( 190 -|q) words compiled 

FORTRAN Input/ 


variable 

Dimension 

Type 

output Description 

A 

1 

DP 

0 

Semimajor axis of the Sun’s 
orbit 

AEROOT 

1 

DP 

0 

Argument of perigee of the Sun 

B1 

2 

DP 

°] 

j Coefficients to transform the 

[ position of the Sun from the 

B2 

2 

DP 

°j 

r • orbital plane to the mean-of- 

epoch equatorial reference 

B3 

2 

DP 

0 J 

' system. 

DAY 

1 

' DP 

I 

Converts days into hours, 
minutes or seconds 

DTOKM 

1 

DP 

I 

Converts distance into 
kilometers 

E 

1 

DP 

0 

Eccentricity of the Sun’s 
orbit 

RAD 

1 

DP 

I 

180 /it 

TWOPI 

1 

DP 

I 

21T 

XJDATE 

• 1 

DP 

I 

Julian day number of the 
desired epoch 

XLSO 

1 

DP 

0 

Mean anomaly of the Sun 

XNS 

1 

DP 

0 

Mean motion of the Sun 
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Figure 29.- PREPS flow chart. 


Page 1 of 1 


100 




77FM50 


3.3.24 PREPT (Subroutine) 

Purpose : Initialize the geopotential coefficients for the Earth (ref. 13) 

Calling sequence : CALL PREPT 

Called by : CONST 

Subroutines/funotions used : None 

Named COMMON: /GEO / RE,CJ2,CS( 187) ,SS( 187) ,IGEOFL 

/RPOOL / FACO,FAC(36) 

Program data : Size = 1031s ( 537 i q ) words compiled 


FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

CS 

187 

DP 

0 

C coefficients of the 
geopotential model in 
the unnormalized form 

IGEOFL 

1 

I 

I/O 

Flag to determine if the 
C and S arrays are set 





= 0 no 
= 1 yes 

SS 

187 

DP 

0 

S coefficients of the 
geopotential model in 


the unnormalized form 
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Figure 30.- PREPT flow chart. 
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3-3-25 PSANS ( Subroutine ) 


Purpose : Analytical theory of the ASOP program. 

Calling sequence : CALL PSANS (NN) 


Called by : ASOP, TIMEPS 

Subroutines/functions used : • DRAG, LON GPP 


Named COMMON : /CBASIC/ 

/CBODY / 
/CONSTW/ 
/DENS / 
/HAMDS / 
/PERTRB/ 
/PS / 
/PSANSV/ 
/PSANS1/ 
/PSANS2/ 
/PSANS3/ 
/PSTIME/ 
/S1STAD/ 
/S1STAV/ 


PI , TWOPI , DEG, R AD , DAY , DTOKM 
XMU , XMUI , SQTMU , SQTMU1 , EPS 
TW03,BY3,BY6, CN- 
BC 10) ,DS1 ,DS2,DC1 , DC2 , 00 , 0C2 , 0S2 
DSF,DSB 
IDRAG, ILONG 
SIG(4),RHO(4) ,TAU 
FACTOR (4 ) , SIGINI ( 8 ) 

ETA1 ,ZETA1 , TOOL, IQL, FAR 

S0M2 , SUM3 , DIFF2 , DIFF3 ,G,H, PSSQRT , PS , QS 

QC , EROOT , X3R00T 

CLO , FAKTPS , TOL 

GC(8) ,P(8) ,Q(8) ,HC(8),QCV(8) 

GIN , HOG , GPH ,BS,FS, GINSQ 


Equivalences : (LS,SIG(4)), (L,RH0(4)), (PHI,RH0(1)) , (DSF( 1) ,W( 1 ) ) , 

(DSB( 1 ) ,W(5) ) , (SI (1 ) ,GAM3( 1 ) ) , (Y(1 ) ,GAM2( 1 ) ) , <Q( 1 ) ,GAM< 1 ) ) , 
(HC(1),DEL3(1)), (GC( 1 ) ,DEL2( 1 ) ) , (P( 1 ) ,DEL( 1 )) 

Program data : Size = 1712g (970 10 ) words compiled 


A description of the mathematical symbols used and their relationship to one an- 
other is given in reference 17; a brief description can also be found in 
Appendixes E and F. Therefore, only a brief mathematical description will be 
given . 


FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

BS 

1 

DP 

0 

b = 1 - H/G 

BY3 

1 

DP 

0 

1/3 

BY6 

1 

DP 

0 

. 1/6 

CN 

1 

DP 

0 

± depending on value i 

DIFF2 

. 1 

DP 

0 

P 2 2 - ° 2 2 

DIFF3 

1 

DP 

0 

p 3 2 . ct 3 2 
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FORTRAN 

variable 

Dimension 

SEe 

Input/ 

output 

Description 

DSB 

4 

DP 

0 

3b/3 k k = 1, 2, 3, 4 

DSF 

4 

DP 

0 

3f/3B k k = 1, 2, 3, 4 

EPS 

1 

DP 

I 

e - 3/2 (V J 2 R e 2 ) 

ETA1 

1 

DP 

0 

sin CT-j 

FACTOR 

4 

DP 

I/O 

Derivatives of the DS 
Hamiltonian and its combina- 
tions (A-j , A 2 , Ag and Ajj in 
appendix F) 

FAK 

1 

DP 

0 

(2L)“ 3/2 

FAKTPS 

1 

DP 

0 

Derivative of the DS 
Hamiltonian (Au in appendix 
F) 

FS 

1 

DP 

0 

f = 1/pq 

G 

1 

DP 

0 

G 

GC 

8 

DP 

0 

3G/3CT k , 3G/8pi_ 
k = 1, 2, 3, 4 

GIN 

1 

DP 

0 

G" 1 

GINSQ 

1 

DP 

0 

G- 2 

GPH 

1 

DP 

0 

G + H 

H 

' 1 

DP 

0 

' H 

HC ' 

8 

DP 

0 

3H/3a k , 3H/3p^ 
k = 1, 2, 3, 4 

HOG 

1 

DP 

0 

H/G 

IDRAG 

1 

I 

I 

Flag to determine if drag cal 
culations are to be included 





= 0 no 
= 1 yes • 

IQL 

1 

DP 

0 

y/V37 

L 

1 

DP 

I/O - 

L = p 4 = a 8 
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FORTRAN 



Input/ 

variable 

Dimension 

Type 

DP 

output Description 

LS 

1 

I/O 

1 = 

NN 

1 

I 

I 

Flag determining if initial- 
izing or computing 





= 0 initializing 





= 1 computing 

00 

1 

DP 

°) 

j Coefficients of Fourier 

series describing density 

0C2 

1 

DP 

° I 

> variation due to J 2 changes- 

1 in height (ref. 7) 

0S2 

1 

DP 

0 J 

1 

P 

8 

DP 

0 

3p/3cr k , 3p/3p k 
k • 1j 2, 3y 4 

PHI 

1 

DP 

I/O 

4> = P-| = a 5 

PS 

1 

DP 

0 

P 

PSSQRT 

1 

DP 

0 

Vp" 

Q 

8 

DP 

0 

3q/3a k , 3q/3p k 
k = 1, 2, 3, 4 

QC 

,1 

DP 

0 

Q 

QCV 

8 

DP 

0 

3Q/3a k , 3Q/3p k 
k = 1, 2, 3, 4 

QS 

1 

DP 

0 

q 

RHO 

4 

DP 

I/O 

p^ , . . .,pl} (see SIG) 

SIG 

4 

DP 

I/O 

a i , . . - , CTj^ ; note the location, 
of SIG and RHO in COMMON. 

This location makes the equiv- 
alence 

Pi = °5> * • • >P4 = a 8 

SIGINI 

8 

DP 

0 

The initial values of the cr’s 

\ 




and p's 

SQTMU 

1 

DP 

I' 

vr 

SQTMUI 

1 

DP 

I 

1/VTT 

SUM2 

1 

DP 

0 

1/2 (Cf 2 2 + P 2 2 ) 
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FORTRAN 

variable 

Dimension 

Type 

SUM3 

1 

DP 

TAU 

1 

DP 

TWOL 

1 

DP 

TW03 

1 

DP 

XMU 

1 

DP 

XMUI 

1 

DP 

ZET1 

1 

DP 


Input/ 

output 

Description 

0 

1/2 (a 3 2 + p 3 2 ) 

. I 

PS Independent variable 

0 

2L . 

0 

2/3 

I 

V 

I 

u- 1 - 

0 

cos cr-j 
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Figure- 31.- PSANS flow charts. 
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Figure 31.- Continued. 
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Figure 31.- Continued. 
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Figure 31.-* Concluded. 
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3.3*26 PST OX (Subroutine) 

Purpose : Transform the PS (Poincare-Similar) elements into the Cartesian coordi 

nates (X,V); the subroutine will also compute the physical time for 
the time iteration stopping procedure. 


Calling sequence : CALL PSTOX (ITIME) 

Called by : ASOP, PREPD, TIMEPS 

Subroutines /functions used: GEOPOT 


Named COMMON: /CARTC / 

/CBODY / 
/PS / 
/PSANS1 / 
/PSANS2/ 
/PSANS3/ 
/PSTIME/ 
/RETRO / 
/RPOOL / 


XI, X2,X3, VI ,V2,V3, TIME, ENERGY, R,RI 

XMU , XMUI , SQTMU , SQTMOI ,EPS 

SIG(4),RH0(i|),TAU 

SSIG1 , CSIG1 , TWOL,XIQL ,FAK 

SUM2 , SUM3 , DIFF2 , DIFF3 , GC ,HC , PSSQRT , PS , QS 

QC,ER00T,X3R00T 

CLO , FAKPS , TOL 

IRO 

ECOSPH , ESINPH , ROP , EMINPH , GCIN , RCAP , RDOT , RCAPDT , XXX ( 6 ) 


Program data : Size = 3318 (217io) words compiled 


FORTRAN 

variable 

Dimension 

Type, 

Input/ 

output 

Description 

CS1G1 

1 

DP 

I 

cos <?i 

ENERGY 

1 

DP 

0 

Total energy of the system 

EROOT 

.1 

DP 

I 

"V 2Lp/p 

FAK 

1 

DP 

I 

(2L)-3/2 

GC 

1 

DP 

I 

Pi - 1/2 (<r 2 2 + P 2 2 ) = G 

IRO 

1 

DP 

I 

Flag to determine if the 


orbit is retrograde 

•= -1 yes 
= 1 no 
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FORTRAN 

variable 

Dimension 

Type 

Input/ 

output Description 

ITIME 

1 

I 

I 

Flag to determine which terms 

are to be computed 

=0 Compute all Cartesian co- 
ordinate elements 

=1 Compute only physical 
time 

=2 Compute only position 
vector 

Note: ITIME is assumed to be 

0 if it is not 1 or 2 

POT 

1 

DP 

I 

Magnitude of the Earth's 
gravitational potential 

PS 

1 

DP 

I 

p (see section 4.2) 

QC' 

1 

DP 

I 

' Q (see section 4.2) 

QS 

1 

DP 

I 

q (see section 4.2) 

R 

1 

DP 

0 

Magnitude of the position 
vector of the satellite 

RHO 

4 

DP 

I 

p-j,...,P4 

RI 

1 

DP 

0 

Inverse magnitude of the posi- 
tion vector of the satellite 

SIG 

4 

DP 

I 

CT 1 1 • • • > ct 4. 

SSIG1 ‘ 

1 

DP 

I 

sin cr-j 

SOM3 

1 

DP 

I 

a 3 2 + p 3 2 

TIME 

1 

DP 

0 

t 

TWOL 

1 

DP 

I 

2 P4 

VI 

1 

DP 

°) 

| Components of the velocity 

[ vector 

V2 

1 

DP 

“I 

^ VI = X component 

i V2 = Y component 

V3 

1 

DP 

0 J 

’ V3 - Z component 
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FORTRAN 

variable 

Dimension 

Type 

Input/ 

output Description 

XI . 

1 

DP 

° 1 

j Components of the position 

vector 

X2 

1 

DP 

° | 

? XI = X component 

X2 = Y component 

X3 

1 

DP 

0 j 

' X3 = Z component 

X3R00T 

1 

DP 

I 

(V^G - Cfg 2 _ p 3 2 )/g 

XMU 

1 

DP 

I 

V 

XMUI 

I 

DP 

I 

U' 1 
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Figure 32.- PSTOX flow chart. 
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3-3.27 RECUR (Subroutine) 

Purpose : Compute, recursively, the sine and cosine of multiples of an angle 

Calling sequence : CALL RECUR (COS, SIN, MAX) 

Called by : POTEXP 

Subroutines/functions used: None 


Named COMMON: None 


Program data : Size = l63g (115 -jq) words compiled 


FORTRAN 

variable 

— “ U 

Dimension 

Type 

Input/ 

output 


Description 

COS 

MAX 

DP 

I/O 

C0S(1) 

cos 

= cosine of initial 
angle 

= The array cos(n6) , 
n=1 ,MAX 

MAX 

1 

I 

I 

Number of terms to be 
generated 

SIN 

MAX 

DP 

I/O 

SIN( 1 ) 
SIN 

= sine of initial 
angle 

= The array sin(n6) , 
n=1 ,MAX 


IIP. 



^ BEGIN ^ 


COMPUTE, 
RECURSIVELY, 
THE SINE AND 
COSINE OF 
MULTIPLES OF 
AN ANGLE 


^ RETURN ^ 


Figure 33.- RECUR flow chart. 
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3.3.28 SACT (Subroutine) 

Purpose : Determine the solar activity coefficients for a given date (ref. 19) 

Calling sequence : CALL SACT (DAY, F10BAR, AKP) 

Called by : INITAL 

Subroutines/functions used : None 

Named COMMON : None 

Program data : Size = 5458 (357io) words compiled 


FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

AKP 

1 

DP . 

0 

Averaged value for the geomag- 
netic index 

DAY 

1 

DP 

I 

Julian day number for which the 
solar activity coefficients are 
desired 

F10BAR 

1 

DP 

0 

Averaged value for the solar flux 
coefficient, Fiq.7 
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Figure 34.- SACT flow chart. 
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3.3.29 SUN (Subroutine) 

Purpose ; Compute the Sun's position analytically (ref. 18 ) 

Calling sequence : CALL SUN (TIME) 

Called by : INITAL 

Subroutines /functions used : MTOECC 

Named COMMON : /CBASIC/ PI,TWOPI,DEG,RAD,DAY,DTOKM 

/SUNPAR/ XNS , XLSO , A , E , AEROOT , B 1 .( 2 ) ,B2(2) ,B3(2) 
/SUNPOS/ XS , YS , ZS , RS , RAS , DECS 

Program data : Size = 1538 (107 -jq) words compiled 


FORTRAN 

variable 

Dimension 

Type 

Input/ 

output Description 

A 

i 

DP 

I 

Semimajor axis of the Sun's orbit 

AEROOT 

i 

DP 

I 

Argument of perigee of Sun 

B1 

2 

DP 

"I 

| Coefficients to transform the 

f position of Sun from orbital 

B2 

2 

DP 

x | 

> plane to mean-of -epoch 

equatorial reference system 

B3 

2 

DP 

I J 


COSEA 

1 

DP 

I 

Cosine of eccentric anomaly 

DECS 

1 

DP 

0 

Declination of Sun 

E 

1 

DP 

I 

Eccentricity of Sun 

EA 

1 

DP 

I 

Eccentric anomaly (rad) 

RAS 

1 

DP 

0 

Right ascension of Sun 

RS 

1 

DP 

0 

Magnitude of the position 
vector of Sun 

SINEA 

1 

DP 

I 

Sine of eccentric anomaly 

TIME 

1 

DP 

I 

Elapsed time of epoch 

TWOPI . 

1 

DP 

I 

2it 

XLSO 

1 

DP 

I 

Mean anomaly of Sun 

XM 

1 

DP 

0 

Mean anomaly (rad) 
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FORTRAN Input/ 


variable 

Dimension 

Type 

output 

DESCRIPTION 

XNS 

1 

DP 

I 

Mean motion of Sun 

XS 

1 

DP 

p ! 

Components of the Sun's position 
vector in the Earth's .inertial 

YS 

1 

DP 

0 

equatorial system 

ZS 

1 

DP 

0 } 
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( BEGIN ^ 


CALL MTOECC 


ANALYTICALLY 
DETERMINE THE 
SUN’S POSITION 
ITS MAGNITUDE, 
THE RIGHT 
ASCENSION OF 
THE SUN AND 
THE DECLINATION 
OF THE SUN 


l 

RETURN ^ 
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Figure 35.- SUN flow chart. 
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3.3.30 TABLE (Subroutine) 

Purpose, ; Generate the table of coefficients for the sine and cosine of 
rP i for the drag disturbing function (refs. 6 and 7) 

Calling sefluence, :, CALL TABLE 


Called by : PREFD 

Subroutines/functions used : None 

Named COMMON: /DBETAS/ B1,B2,B3,B4 

/DENS / COEFF(IO) ,DSi,DS2,DC1 ,DC2,HM,HC2,HS2,HC1 ,HS1 
/DTABLE/ TT( 12) 

/FORSAV/ Z(6),ZC(9,6),ZS(9,6),DZ(6),DZC(9,6),DZS(9,6) 
/PS / SIG(4) , RH0(4) ,TAU 


Program data 

: Size = 

I423g (787 

FORTRAN 

variable 

Dimension 

Type 

B1 

1 

DP 

COEFF 

10 

DP 

DC1 

1 

DP 

DC2 

1 

DP 

DS1 

1 

DP 

DS2 

1 

DP 

DZ 

6 

DP 

DZC 

(9,6) 

DP 

DZS 

(9,6) 

DP 

HC1 

1 

DP 

HC2 

1 

DP 

HS1 

1 

DP 

HS2 

1 

DP 

RHO 

4 

DP 


Input/ 

output 

I 

I 

I 
I 
I 
I 
0 
0 
0 
I 
I 
I 
I 
I 


Description 


a 


Coefficients of the density 
model (ref. 7) 


Coefficients of the Fourier 
series describing the diurnal 
bulge 


Coefficients of the Fourier 
series expansion of the 
disturbing function in 
PS elements (ref. 7) 


Coefficients of the Fourier series 
describing the density variation 
due to J2 changes in height 


PS elements P-j, (>£, P 3? P^ 
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FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

SIG 

4 

DP 

I 

PS elements <*i , o 2 , <Jg , aij 

TT 

12 

DP 

0 

Table of the averaged Fourier 
series coefficients 

Z 

6 

DP 

0 

Xq 1 Coefficients of the Fourier 
. 1 series expansion of the 

zc 

(9,6) 

DP 

0 

> disturbing function in 
I PS elements (ref. 7) 

ZS 

(9,6) 

DP 

0 
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c 


BEGIN 


GENERATE TABLE 
OF COEFFICIENTS 
OF no-j FOR THE 

DRAG DISTURBING 
FUNCTION 


RETURN 


3 


Figure 36.- TABLE flow chart. 
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3.3.31 TIMEPS (Subroutine) 

Purpose : Iteration procedure to stop the PS elements at a desired value of 

the physical time 

Calling sequence : CALL TIMEPS (TFIN,ISET) 


Called by : ASOP 

Subroutines/functions used : PSANS , PSTOX 


Mamed COMMON : /CARTC / 

/ CBASIC/ 
/PS / 
/PSANS 1/ 
/PSANS2/ 
./PSTIME/ 


X( 6 ) , TIMEP , ENERGY , R , RI 
PI , TWOPI , DEG , RAD , DAY , DTOKM 
SIG ( 8 ) , -TAU , TAUMAX 
SSIG1 , CSIG1 

SUM2 , SUM3 , DIFF2 , DIFF3 , G , H ,PSSQRT ,PS , QS 
CLO,FAKT,TOL 


Program data: 


Size = (292-jg) words compiled 


FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

CLO 

1 

DP 

I 

Initial value of 

CSIG1 . 

1 

DP 

0 

cos 

DAY 

1 

DP 

I 

Value to convert days into 
hours, minutes, or seconds 

FAKT 

1 

DP 

I 

FAKTPS from PSANS = A 4 (see 
appendix F) 

I SET 

1 

DP 

0 

Flag to determine - if the 
force model must be updated 
= 0 no (Force model still valid) 
= 1 yes (Outside- range of 

validity for force model 
update state vector to 
Tma x and reinitialize 
force model 

QS 

1 

DP 

I 

q (sec. 4.2) 

RI 

1 

DP 

I 

Inverse of the magnitude of the 
position vector of the satellite 

SSIG1 ' 

1 

DP 

I- 

sin CT-| 
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FORTRAN 

variable 

Dimension 

. Type 

Input/ 

output 

Description 

TAU 

1 

DP 

o - 

New value of the independent 
variable ( T ) 

TAUMAX 

1 

DP 

I 

Range of validity for the force 
,-model 

TFIN 

1 

DP 

i' 

Final time desired for stopping 
the iteration 

TIMEP 

1 

DP 

I 

* 

Computed value of the physical time 

TOL 

1 

DP 

0 

Allowable TOLerance between TFIN 
and TIME that will stop the 
iteration 





|TFIN - time! < TOL 

TWOPI 

1 

DP 

I 

• 21T 
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3.3-32 TIMEXP (Subroutine) 

Purpose : Compute coefficients of the expansion of the time equation 

Calling sequence r CALL TIMEXP (XNU,ECC2,XMCAP,NMAX) 

Called by : POTEXP 

Sub routines /functions used : None 

Named COMMO N: None 

Program data : Size = 200g ( 1 28 1 0 ) words compiled 


FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

ECC2 

i 

DP 

I 

e 2 

NMAX 

i 

I 

I ■ 

Number of terms to compute 

XMCAP 

(7,18) 

DP 

0 

Fourier coefficients of the 
time expansion 

XNU 

1 

DP 

I 

Ratio of frequency of rotation 
of satellite to rotation of 
Earth 
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Figure 38„- TIMEXP flow chart* 
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3.3-33 XTOPS. (Subroutine) 

Purpose : Transform the Cartesian coordinates (X,V) into the PS (Poincare- 

Similar) elements (a,p) 


Calling sequence : CALL XTOPS 

Called by : ASOP, INPUT 

Subroutines/f unctions used: GEOPOT 


Named COMMON : /CARTC / 

/CBODY / 
/PS / 
/PSANS1/ 
/PSANS2/ 
/PSANS3/ 
/PSTIME/ 
/RETRO/ 
/RPOOL/ 


XI , X2 , X3 , VI , V2 ,V3 , TIME , ENERGY , R , RI 

XMU , XMUI , SQTMU , SQTMUI , EPS 

SIG(4),RHO(4),TAU 

SSIG 1 , CSIG 1 , TWOL , XIQL , FAX 

SUM2 , SUM3 , DIFF2 , DIFF3 , GC , HC , PSSQRT , PS , QS 

QCAP 

CLO , FAKTPS , TOL 
IRO 

G 1 SQ , G2SQ , ESINPH , EROOT , EMINPH , SQTGHI , RCAP , RDOT , 
ZCAP1 , ZCAP2 , ECOSPH ,XXX( 4 ) 


Equivalence : (HC,G3), (PHIC,RH0(1)) , (XL,RHO(4)) 

Program data : Size = 576g (382 -|q) words compiled 

FORTRAN Input/ 

variable Dimension Type output Description 

(All equations are described in section 4.1) 


CLO 

1 

DP 

0 

Initial value of cty 

CSIG1 

1 

DP 

0 

cos a-j 

DIFF2 

1 

DP 

0 

p 2 2 - a 2 2 

DIFF3 

1 

DP 

0 

p 3 2 - a 3 2 

ENERGY 

1 

DP 

0 . 

Total orbital energy 

FAK 

1 

DP 

0 

(2L)“ 3/2 

GC 

' 1 

DP 

0 

G = ^G x 2 + G y 2 + G z 2 

HC 

1 

DP 

0 • 

H = G z 

IRO 

1 

I 

0 

Flag to determine if the orbit 
is retrograde. 


= -1 yes 

= 1 no 
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FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

POT 

1 

DP 

I 

PS 

1 

DP 

0 

PSSQRT 

1 

DP 

0 

QCAP 

1 

DP 

0 

QS 

1 

DP 

0 

R 

1 

DP 

0 

RHO 

4 

DP 

0 

RI 

1 

DP 

0 

SIG 

4 

DP 

0 

SSIG1 

1 

DP 

0 

SQTMUI 

1 

DP 

I 

SUM2 

1 

DP 

0 

SUM3 

1 

DP 

0 

TAU 

1 

DP ' 

0 

TIME 

1 

DP 

I 

TWOL - 

1 

DP 

0 

VI 

1 

DP 

1 ) 

V2 

1 

DP 

1 

V3 

1 

DP 

I ) 

XL 

1 

DP 

0 

XI 

1 

DP 

1 ) 

X2 

1 

DP 

1 

X3 

1 

DP 

I ) 


Description 


Magnitude of the Earth’s 
gravitational potential 

P 

VP" 

Q 

Q 

Magnitude of the position vector 
of the satellite (r) 

pi -> p4 in RH0(1) + RH0(4) 

Inverse magnitude of the position 
vector of the satellite ( 1 /r) 

Oi ♦ CT 4 in SIG(1) * SIG(4) 

sin C-] 

1/VjT 

1/2 (a 2 2 + P 2 2 ) 

1/2 (CT 3 2 + P 3 2 ) 

PS independent variable; 
initially set to -zero (rad) 

Initial physical time 


2 L = 2p 4 = 2cfg 


Components of the velocity vector 

VI = X component 
V2 = Y component 
V3 = Z component 

P 4 

Components of the position vector 
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FORTRAN 

variable 

Dimension 

T ype 

Input/ 

output 

Description 

XIQL 

1 

DP 

0 

y/ V2L 

XMU 

1 

DP 

.1 

y-' 
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Figure 39.- XTOPS flow chart. 
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3-3-3^ XYZASI (Subroutine) 

P urpose ; Transform the Cartesian coordinates (X,V) into the Keplerian elements 
(a,e,i,w,fl,M) 


Calling sequence .: CALL XYZAEI 


Called by : OUTPUT 


Sub routines/ functions use d: None 


Named COMMON : /CARTC / 

/CBASIC/ 
/CBODY / 
/KEPLER/ 
/RPOOL / 


XI , X2 , X3 , V 1 , V2 , V3 , TIME , ENERGY , R , RI 
PI , TWOPI , BEG , RAD , DA Y , DTOKM 
XMU , XMUI , SQTMU , SQTMUI , EPS 
A , E , XI , OMEGA , XNODE , XM 

VSQ , RRDOT , ECOSE , ESINE , P , H , WX , WY , WZ , EA , RCOSF , RSINF , 
COSWF ,SINWF , TEMP, RCOSL ,RSINL , XXX{ 7 ) 


Program data : Size = *l06g (262 -jq) words compiled 

Considers only elliptic motion 


FORTRAN 

variable 

Dimensi on 

Type 

Input/ 

output Description 

A 

i 

DP 

0 

Semimajor axis of the orbit 

E 

i 

DP 

0 

Orbit eccentricity (e) 

OMEGA 

i 

DP 

0 

Argument of pericenter (co) 

R 

i 

DP 

0 

Magnitude of the position vector 
of the satellite 

SQTMUI 

i 

DP 

I 

1/ VTT 

VI 


DP 

I N 

i 

v Components of the velocity vector 

V2 

i 

DP 

1 

I 

1 

/ VI = X component 

1 V2 = Y component 

V3 

i 

DP 

I , 

' V3 = Z component 

XI 

i 

DP 

' ) 

, Components of the position vector 

X2 

i 

DP 

' 

' XI = X component 

X2 = Y component 

X3 ' 

i 

DP 

I ) 

X3 = Z component 

XM 

i 

DP 

0 

Mean anomaly (rad) 


134 



77FM50 


FORTRAN 

variable 

Dimension 

Type 

Input/ 

output 

Description 

XMU 

i 

DP 

I 

Central body gravitational 
constant (y) 

XNODE 

i 

DP 

0 

Argument of the ascending 
node (ft) (rad) 
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Figure 40.- XYZAEI flow chart. 
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3.4 LABELED COMMON: 

Notation: R - Real variable 

I - Integer variable 

S - Single, .pr.ecision- 
D - Double precision 


Type 


Precision 


/CARTC/; Cartesian coordinates of the satellite's position 


Location 

Name 

Dimension 

Type/ 

precision 

Description 

1-3 

X 

3 

R/D 

Cartesian coordinates for the 
position of the satellite; 
X,Y,Z 

4-6 

V 

3 

R/D 

Velocity vector of the satel- 
litej V x ,V y ,V z 

7 

TIME 

1 

R/D 

Elapsed time (hr, min, or sec) 

8 

ENERGY 

1 

R/D 

Total orbital energy 

9 

R 

1 

R/D 

Magnitude of the position 
vector of the satellite 

10 

RI 

1 

R/D 

Inverse magnitude of the posi- 
tion vector of the satellite 


In subroutines : MAIN, AEIXYZ, ASOP, DENSTY, GEOPOT, INPUT, 

OUTPUT, PREPD, PSTOX, TIMEPS, XTOPS, XYZAEI 


/CBASIC/; 

Conversion 

constants 




Location 

Name 

Dimension 

Type/ 

precision 


Description 

1 

PI 

1 

R/D 

IT 


2 

TWOPI 

1 

R/D 

2-jt 


3 

DEG 

1 

R/D 

it/ 180 


4 

RAD 

1 

R/D 

180/tt 
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Location 

Name 

Dimension 

Type/ 

precision 

Description 

5 

DAY 

i 

R/D 

Converts days into hours, 
minutes, or seconds 

6 

DTOKM 

i 

R/D 

Converts distance into 
kilometers 

In subroutines: MAIN, ASOP, CONST, DENSTY, GEOPOT, INITAL , INPUT, OUTPUT, 

POTEXP, PREPD, PREPS, PSANS, SUN, TIMEPS, XYZAEI 

/CBODY/; 

Gravitational 

variables 



Location 

Name 

Dimension 

Tuype/ 

precision 

Description 

1 

XMU 

1 

R/D 

Gravitational constant of 
the central body (y) 

2 

XMUI 

1 

R/D 

i/u 

3 

SQTMU 

1 

R/D 

VTT 

4 

SQTMUI 

1 

R/D 

VT7 ~v~ 

5 

EPS 

1 

R/D 

e = 3/2 (y J 2 Re^) where 


J 2 = ^2 geopotential 
coefficient 


R g = Equatorial radius of 
the central body 

In subroutines : AEIXYZ, CANFOR, CONST, GEOPOT, LONGPP, POTEXP, 

PREPD, PSANS, PSTOX, XTOPS, XYZAEI 
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/CONSTW/ 





Type/ 


Location 

Name 

Dimension 

precision Description 

i 

TN03 

1 

R/D. 

2/3 

2 

BY 3 

1 

R/D 

1/3 

3 

BY6 

1 

R/D 

1/6 

4 

CN 

1 

R/D 

+1 depending on NN 
= +1 initializing (NN=0) 





= -1 computing (NN/0) 

In subroutines: LONGPP 

, PSANS 



/CPRINT/ 



Type/ 


Location 

Name 

Dimension 

precision Description 

1 

PRINT 

1 

R/D > 

j 

2 

IPRINT 

1 

1 1 

1 Input parameters, see section 

> 2.1.1; Note: IPRINT is set 

3 

IPSPRT 

1 

1 1 

- to 3 if re-entry condition 

i exists 

4 ' 

IUNIT-S 

1 

1 J 

' 

In subroutines : MAIN , 

CONST, INPUT, 

OUTPUT 


/D ATMOS/; 

Atmospheric parameters 






Type/ 


Location 

Name 

Dimension 

precision Description 

1 

FBAR 

i 

R/D 

Averaged value for_the solar 
flux coefficient, F-|o y 

2 

XKP 

i 

R/D 

Averaged value for the geo- 
magnetic index, Kp 

3 

SLDAY 

i 

R/D 

Magnitude of the seasonal 
latitudinal density variation 

4 

SADAY 

i 

R/D 

Magnitude of the semiannual 
density variation 

In subroutines : DENSTY 

, INITAL 
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/DATM01/; Atmospheric parameters 


Location 

Name 

Dimension 

Type/ 

precision 

Description 

1 

SRAB 

1 

R/D 

Sine of bulge right ascension 

2 

CRAB 

1 

R/D 

Cosine of bulge right 
ascension 

' 3 

SDEC 

1 

R/D 

Sine of bulge declination 

4 

CDEC 

1 

R/D 

Cosine of bulge declination 

5 

RR 

1 

R/D 

Magnitude of the diurnal change 
in the exospheric temperature 

6 

TC 

1 

R/D 

Nighttime minimum of the global 
exospheric temperature (°K) 

7 

TG 

1 

R/D 

Variation in the exospheric 
temperature due to geomagnetic 
activity (°K) 

In subroutines: DENSTY 

, INITAL 



/DBETAS/; 

Powers of beta 



Location 

Name 

Dimension 

Type/ 

precision 

Description 

1 

BETA1 

1 

R/D 

$ 

2 

BETA2 

1 

R/D 

&2 

3 

BETA3 

1 

R/D 

&3 

4 

BETA4 

1 

R/D 


In subroutines : CANFOR 

, PREPD, TABLE 



/DCOEFF/; 

Coefficients 

for the Jacchia 71/Lineberry atmospheric density model 

Location 

Name 

Dimension 

Type/ 

precision 

Description 

1-27 

A 

(3,3,3) 

R/D 

Parameters for determining the 
base altitude 
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Type/ 


Location 

Name Dimension 

precision 

Description 

28-54 

B 

(3,9) 

R/D 

Parameters for determining 
the Too density profile 

55-81 

C 

(3,9) 

R/D 

Parameters for computing 
the annual variation 

82-93 

D 

(3,4) 

R/D 

Parameters for computing the 
seasonal latitudinal variation 

In subroutines: DENSTY, 

INITAL 



/DENS/; 

ASOP atmospheric 

density model parameters 



Location 

Name 

Dimension 

Type/ 

precision Description 

1-10 

B 

10 

R/D 

Coefficients of the atmo- 
spheric density model a 

11 

* DS1 

1 

R/D ^ 


12 

DS2 

1 

R/D | 

^ Coefficients of the Fourier 

> series describing the diurnal 

13 

DC1 

1 

R/D I 

^ bulge a 

14 

DC2 

1 

R/D J 


15 

00 

1 

R/D > 

Coefficients of the Fourier 

16 

0C2 

1 

R/D | 

series describing the density 
variation due to short period 

17 

0S2 

1 

R/D 

> J 2 changes in height 3 

18 . 

0C1 

1 

R/D 1 

* 

19 

0S1 

1 

R/D , 


In subroutines : LONGPP 

, PREPD, PSANS 

, TABLE 
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/DETE/; First-order long period parameters 


Location 

Name 

Dimension 

Type/ 

precision Description 

i 

SHAT 

1 

R/D 

First-order, long-period 





generating function b 

2 

SHATP 

1 

R/D 1 

j Derivatives of the first-order, 

3 

SHATE2 

1 

R/D 

i long-period generating function 

4 

SHATB 

1 

R/D 

> with respect to p,e^,b,xi and 

5 

SHATXI 

1 

R/D | 

ih b 

6 

SHATPI 

1 

R/D , 


In subroutines: DETERM, LONGPP 



/DRAG/; Input parameters 






Type/ 


Location 

Name 

Dimension 

precision Description 


1 CD 1 

2 AREA 1 

3 XMASS 1 

4 CDRAG 1 



Input parameters, see section 
2 . 1.1 


In subroutines : CONST, INPUT, PREPD 


/DRAG1/; Drag functions 


Location 


1-8 


9 


Name 

Dimension 

Type/ 

precision 

Description 

CFORCE 

8 

R/D 

Drag force defined in PS 
elements 

T4 

1 

R/D 

Magnitude of the quadratic 
variation in the mean anomaly 


^Reference 7, pp. 18-22. 
^Reference 8. 
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Type/ 


Location 

Name 

Dimension 

precision 

Description 

10 

TLINER 

1 

R/D 

Linear change in time due to 
drag 

11 

TEMPO 

1 

R/D 

Second order correction for 
density due to drag 

In subroutines : DRAG, 

PREPD 



/DTABLE/; Averaged drag functions 






Type/ 


Location 

Name 

Dimension 

precision 

Description 

1-12 

TT 

12 

R/D 

Table of averaged Fourier 
series in 

In subroutines : CANFOR 

, TABLE 



/ECC/; Eccentricity parameters 


- 




Type/ 

• 

Location 

Name 

Dimension 

precision 

Description 

1 

ES 

1 

R/D 

e 

2 

ESSQ 

1 

R/D 

e 2 

In subroutines: DETERM 

, FPRIME, LON GPP 


/END/; Input 

stopping parameters 






Type/ 


Location 

Name 

Dimension 

precision 

Description 

1 

STOP 

1 

R/D 1 

Input parameters, see section 
2.1.1; Note: ISTOP is reset 

2 

ISTOP 

1 

I } 

to 3 if reentry condition 
exists 


In subroutines : MAIN, CONST, INPUT, OUTPUT, PREPD 
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/EPOCH/ ; Input dating parameters 


Location 

Name 

Dimension 

Type/ 

precision 

. Description 

1-6 

DATE 

6 


R/D 

Input parameter, see section 
2.1.1 

7 

XJDATE 

1 


R/D 

Julian date corresponding to 
DATE 

In subroutines: CONST, 

INITAL,, 

INPUT, OUTPUT, PREPS 

/EXPCOF/; ; 
and sine 

Binomial coefficients 

and 

Fourier coefficients of the powers of cosine 

Location 

Name 

Dimension 

Type/ 

precision 

Description 

1-200 

A 

200 


R/D 

Array containing the binomial 
coefficients 

201-400 

B 

200 


R/D 

Array containing the Fourier 
coefficients for cosine and 
sine raised to a power 

401 

NDEX0 

1 

- 

I 

Zero Index to the NDEX array 

402-419 

NDEX 

18 


I 

Array of pointers to the A and 
B coefficients 

420 

IEXPFL 

1 


I 

Flag to determine if the A and 
B arrays have been computed 


= 0 no 
= 1 yes 

In subroutines : COEFF, DETERM, FPRIME, POTEXP 
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/FORSAV/; Fourier series expansion coefficients 


Location 

Name 

Dimension 

Type/ 

precision 

Description 

•1-6 

z 

6 

R/D 

x 1 N 

o 


7-60 

zc 

(9,6) 

R/D 

X 1 i 







Coefficients of the 

61-114 

zs 

(9,6) 

R/D 

^ I 

[ Fourier series expansion 





J \ 

A 

\ of the disturbing 

115-120 

\ 

DZ 

6 

R/D 

X 0 1 

' function in PS elements a 

121-174 

DZC 

(9,6) 

R/D 

< X 


175-258 

DZS 

(9,6) 

R/D 

A 

1 


In subroutines: PREPD, TABLE 


/ 


/FP/ ; Second-order 

zonal parameters 






Type/ 



Location 

Name 

Dimension 

precision 

Description 

i 

FHAT 

1 

R/D 

Second-order zonal Hamiltonian* 5 

2 

FHATP 

1 

R/D 'I 

| Derivatives of the second-order 




| 

f zonal Hamiltonian with respect 

3 

FHATE 2 

1 

R/D , 

r to p 

, e 2 and b ° 

4 

FHATB 

1 

R/D J 

1 



In subroutines : FPRIME, LONGPP 


/GEO/; Geopotential coefficients for the Earth 


Location 

Name 

Dimension 

Type/ 

precision 

Description 

1 

RE 

i 

R/D 

Equatorial radius of the 
central body (Earth) 

2 

CJ2 

i 

R/D 

J 2 geopotential coefficient 
of the central body (Earth) 


^Reference 7, pp. 43-45. 
Reference 8. 
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Location 

Name 

Dimension 

Type/ 

precision 

Description 

3-189 

CS 

187 

R/D 

C coefficients of the 
geopotential model in the 
unnormalized form 

190-376 

SS 

187 

R/D 

S coefficients of the 
geopotential model in the 
unnormalized form 

377 

IGEOFL 

1 

R/D 

Flag to determine if C and 
S arrays are set 

= 0 no 
= 1 .yes 


In subroutines : ASOP, CONST, DETERM, FPRIME , GEOPOT, POTEXP, PREPD , PREPT 


/GMTROT/; Greenwich Meridian rotational parameters referenced to- a desired epoch 


Location 

Name 

Dimension 

Type/ 

precision 

Description 

1 

WE 

1 

R/D 

Rotational rate of the Earth 

2 

THETAO 

1 

R/D 

Initial hour angle of the Earth 

In subroutines: CANFOR, CONST, GEOPOT. POTEXP. 

PREPD 

/HAMDS/ ; 

Hamiltonian 

derivatives 



Location 

Name 

Dimension 

Type/ 

precision 

Description 

1-4 

DSF 

4 

R/D 

3f/38 K k=1, 2, 3, 4 

5-8 

DSB 

4 

R/D 

3b/9B K - k=1, 2, 3, 4 


In subroutines : LONGPP , PSANS 
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/KEPLER/; Keplerian elements of the satellite's orbit 

Type/ 


Location 

Name 

Dimension 

precision 

Description 

i 

A 

1 

R/D 

Semimajor axis of the orbit 

2 

E 

1 

R/D 

(a) 

Eccentricity (e) 

3 

XI 

1 

R/D 

Orbital inclination to the 

4 

OMEGA 

1 

R/D 

Equator (i) 

Argument of pericenter ( u ) 

5 

XNODE 

1 

R/D 

Argument of the ascending 

6 

XM 

1 

R/D 

(f2) node 

Mean anomaly (M) 

In subroutines: AEIXYZ, INPUT, OUTPUT, XYZAEI 


/PERTRB/; 

Location 

Input perturbation flags 
Name Dimension 

Type/ 

precision 

Description 

1 

IDRAG 

1 

I 

Input parameters, see section 

2 

ILONG 

1 

I 

2.1.1 

In subroutines: ASOP, 

CONST, GE0P0T, 

INPUT, PREPD , PSANS 

/PS/; PS 

elements and 

independent variable 


Location 

Name 

Dimension 

Type/ 

precision 

Description 

1-4 

SIG 

4 

R/D 

PS elements d-], c^, d^, cty 

5-8 

RHO 

4 

R/D 

PS elements p-| , P2, P3, P4 

9 

TAU 

1 

R/D 

Independent variable of 

10 

TAUMAX 

1 

R/D 

the PS elements (t) 

Range of validity for the force 


model; When t exceeds this 
value ( T max ), the force model 
must be. reinitialized 
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Type/ 

Location Name Dimension precision Description 

^ TAUINT 1 R/D Initial value of T for 

which the force model is 
valid; initially set to 0 

In subroutines : MAIN, ASOP, CANFOR, DRAG, LONGPP, OUTPUT, POTEXP, 

PREPD, PSANS, PSTOX, TABLE, TIMEPS, XTOPS 


/PSANSV/; PS formulation variables 

Type/ 


Location 

Name 

Dimension 

precision 

Description 

1-4 

FACTOR 

4 

R/D 

Derivatives of the DS 
Hamiltonian and its combin- 
ations (A 1 , A 2 , Aj and 
Ajj in appendix F) 

5-12 

SIGINI 

8 

R/D 

Initial values of the cr's 
and p’s 

In subroutines: LONGPP, 

PREPD, PSANS 



/PSANS1/; 

PS formulation 

. variables 



Location 

Name 

Dimension 

Type/ 

precision 

Description 

1 

SSIG1 

1 

D/P 

sin cr-j 

2 

CSIG1 

1 

D/P 

cos a*j 

3 

TWOL 

1 

D/P 

2L 

4 

XIQL 

1 

D/P 

U/V2L 

5 

FAR 

1 

D/P 

(2L) _3/2 


In subroutines : CANFOR, LONGPP, POTEXP, PREPD, PSANS, PSTOX, .XTOPS 
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/PSANS2/; PS 

formulation 

variables 



Location 

Name Dimension 

Type/ 

precision 

Description 

1 

SUM2 

1 

D/P 

(02 2 + P2 2 ) /2 

2 

SUM3 

1 

D/P 

(CT3 2 + P3 2 )/2 

3 

DIFF2 

1 

D/P 

P2 2 - Cf2 2 

1} 

DIFF3 

1 

D/P 

P3 2 ~ a 3 2 

5 

G 

1 

D/P 

•G 

6 

H 

1 

D/P 

H 

7 

PSSQRT 

1 

D/P 

Vp 

8 

PS 

1 

D/P 

P 

9 

QS 

1 

D/P 

q 

In subroutines : CANFOR , 

TIMEPS, 

DETERM, FPRIME , LONGPP 
XTOPS 

, POTEXP, PREPD, PSANS, PSTOX, 

/PSANS3/; PS 

formulation 

variables 



Location 

Name 

Dimension 

Type/ 

precision 

Description 

1 

QC 

1 

R/D 

Q 

2 

EROOT 

1 

R/D 

V2L p/u 

3 

X3R00I 

1 

R/D 

(A/HG - 03 2 - P3 2 ) /G 

In subroutines : LONGPP , 

POTEXP, PREPD, PSANS, 

PSTOX, XTOPS 
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/PSTIME/; PS parameters used for stopping on a specific final time- 


Location 

Name 

Dimension 

Type/ 

precision 

Description 

1 

CLO 

i 

R/D 

Initial value of or^ (set 
when initializing, NEWX =0) 

2 

FAKTPS 

i 

R/D 

a<J2j/9T 

3 

TOL 

i 

R/D 

Tolerance criteria for the 
iteration stopping procedure 

In subroutines : MAIN , 

•PSANS, PSTOX, 

TIMEPS, XTOPS 

/RETRO/; 

Retrograde parameter 



Location 

Name 

Dimension 

Type/ 

• precision 

Description 

1 

IRO 

1 

I 

Flag to determine retrograde 


orbit 

= -1 yes, it is a retrograde 
orbit 

= 1 no, it is a retrograde 

orbit 


In subroutines: CANFOR, DETERM, POTEXP, PREPD, 

PSTOX, XTOPS 

/RP00L/,/RP00LA/; Temporary variables 



Location Name Dimension 

Type/ 

precision 

Description 

Variable — 1393 

R/D 

Temporary real variables; 
this COMMON block is used 
to help save storage within 
the ASOP program 


In subroutines : AEIXYZ, DETERM, FPRIME , GEOPOT, LONGPP, POTEXP, PREPT, PSTOX, 

XTOPS , XYZAEI 
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/SUNPAR/; Orbit' parameters of the Sun 


Type/ 


Location 

Name 

Dimension 

precision Description 

i 

XNS 

1 

R/D 

Mean motion of the Sun 

2 

XLSO 

1 

R/D 

Mean anomaly of the Sun 

3 

A 

1 

R/D 

Semimajor axis of the Sun 

4 

E 

1 

R/D 

Eccentricity of the Sun's orbit 

5 

AEROOT 

1 

R/D 

Argument of perigee of the Sun 

6-7 

B1 

2 

R/D N 

j Coefficients to transform the 

position of the Sun from the 

8-9 

B2 

2 

R/D 

> orbital plane to the mean- 

of-epoch equatorial reference 

10-11 

B3 

2 

R/D J 

' system 

In subroutines : PREPD , 

PREPS, SUN 



/SUNPOS/; 

Sun's position referenced 

to an input epoch 




Type/ 


Location 

Name 

Dimension 

precision Description 

1 

XS 

1 

R/D 

Position vector of the Sun 

2 

YS 

1 

R/D 

in Earth's inertial equato- 
rial system (XS, YS, ZS) 

3 

ZS 

1 

R/D 


4 

RS 

1 

-R/D 

Magnitude of the position 
vector of the Sun 

5 

RAS 

1 

R/D 

Right ascension of the Sun 

6 

DECS 

1 

R/D 

Declination of the Sun 


In subroutines : INITAL, PREPD, SUN 
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/S1STAD/ 

Location 

Name 

Dimension 

Type/ 

precision 

Description 

1-8 

GC 

8 

R/D 

3G/3a k , 3G/3p k 

k=1 ,2,3,4 

9-16 

P 

8 

R/D 

3p/3cr k , 3 p/3p k 

k-1 ,2,3,4 

17-24 

Q 

8 

R/D 

3q/3a k ? 3q/3p k 

k=1,2,3,4 

25-32 

HC 

8 

R/D 

3H/aa k , 3H/3p k 

k=1 ,2,3,4 

33-40 

QCV 

8 

R/D 

3Q/3a k , 3Q/3a k 

k=1 ,2,3,4 


In subroutines : LONGPP , PSANS 


/S1STAV/ 

Location 

Name 

Dimension 

Type/ 

precision 

Description 

1 

GIN 

i 

R/D 

G -1 

2 

HOG 

i ' 

R/D 

H/G 

3 

GPH 

i 

R/D 

G+H 

4 

BS 

i 

R/D 

b = 1 - H/G 

Location 

Name 

Dimension 

Type/ 

precision 

Description 

5 

FS 

1 

R/D 

f = 1/pq 

6 

GINSQ 

1 

R/D 

G- 2 

In subroutines: DETERM, 

/TESS/ ; Input parameters 

FPRIME, LONGPP, PSANS 

for the geopotential model (zonal and tesseral terms) 

Location 

Name 

Dimension 

Type/ 

precision 

Description 

/ 

i 

NMAX 

1 

I 

Input parameters; see section 

2 MMAX 

In subroutines: GEOPOT, 

1 I 

INPUT, POTEXP 

2.1.1 
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/XIPSI/ 

Location Name Dimension 

• 1 XII 1 

2 PSI1 


Type/ 

precision Description 

R/D Xl = e s: *- n * G0S S 
I = inclination 
g = argument of perigee 

R/D t|>i = e sin I sin g 
I = inclination 
g = argument of perigee 


In subroutines; DETERM , LONGPP 
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4.0 PS ELEMENT FORMULATION 

The PS (Poincare -Similar) element formulation is described in the following 
subsections . An exact development of this element set is described in reference 
5, and a description of the variables can be found in Appendix E. 

4.1 TRANSFORMATION FROM CARTESIAN COORDINATES TO PS ELEMENTS (XTOPS) 

The transformation into the PS element set (a,p) from a given set of Cartesian 
coordinates (X,V,r) is accomplished within the subroutine XTOPS. It is 
assumed that these coordinates are the initial conditions so that the physical 
time t and the independent variable t both equal zero. The PS elements 
are then computed using the equations 


cr. = tan 


_! ( X x + R cr 3 

fy ♦ a *P3, 


cr„ = 


Z 2 cos cr.j - Z.j sin a. 


» 1 = G --^G 2 + 2r 2 V 


U 


V^P4 


P 2 = Z 2 sin + Z 1 cos a 1 


O = -2G x />/2(G + G z ) 


p = 2Gy/^2(G + G z ) 


CT 4 = t - 


U 


(2p,) 


3/2 


(e - i) 


P 4 


= - - - V 2 -V, where V is the 


r 2 


perturbing potential 


- - Z ? Q 
P * 


V7T7 


The other required relations are: 

^x ” “ ^z^y » ^y = ^z^x ** ^x^z > ®z = ^x^y ~ 

G =^G X 2 + G y 2 + G z 2 

P = 1/p (g - Pl + u/V^PlT) 2 t q = G - 1/2 Pl + ]i/2^2^ 
R # = rR/2G , R = 2GX z /r>/2(G + G z ) 

Z 1 = (p/r - 1)/Q , Z 2 = rp/[Q(2q - G)] 
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Q 

r 

E 


1/U [p 4 (2li/>/2p^ + G z - G)] 1/2 

x*v 

--- Vi. - e 2 .=. p ^2p 4 /U 

<j> = -2 taiT^ZgQ/O +Vl - e 2 + Z,Q)] 


This transformation is performed only once for a given set of Cartesian 
coordinates unless the desired value of X is greater than the maximum value 
for which the force model is still valid . This transformation is repeated 
whenever the force model is reinitialized. 

4.2 TRANSFORMATION FROM PS ELEMENTS TO CARTESIAN COORDINATES (PSTOX) 

This transformation is performed when any intermediate printout is desired 
or when the final condition is met. Therefore, this transformation is coded 
with emphasis on the speed of calculation. 

The Cartesian coordinates defined in terms of the PS elements are given by 
the equations : 


X x = r cos cr^ - R*CTg 

* 

Xy = r sin cr.j - R Pg 

X = R* V2(G + H) 
z 

The other necessary relationships are 

r = p/(1 + e cos (Jj ) , R* = rR/2G ' 

r = e sin <j>/p[2q - p 1 + 1/2 (pg 2 + cr 2 2 )] 

R* = (Rr + Rr)/2G , R = --(p 0 cos cr, - sin a.) 

r 2 i 1 3 t 

2 

R = p^ sin a 1 + cos ct 1 , G = p^ - 1/2 (p 2 + O 


, V = r cos a 1 - - sin a 1 - 


V„ = r sin a, - - cos cr. - 


1 


1 


V = R* "V2(G + S5 

Z 


2 


R*a 


R*p 


3 


3- 
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H = G - 1/2 (p^ + a 3 2) , q = _ 1/2 (p^ + a 2 _ 

p = 1/u [y/>/ 2p 4 - 1/2 (p 2 2 + cr 2 2 )J 

e sin <j) = Q(p 2 sin a ^ + a 2 cos a^) 
e cos- (j> = Q(p 2 cos a 1 - a 2 sin c^) 

Q = 1/p {p 4 [2p/>/2^ - 1/2 (p 2 2 + a /)]} 1/2 
The physical time t is computed using 

t = o n + ~ * “ - e sin A Vi - e 2 J (4.2a) 

^ (2p 4 ) 3/2 \ P / 

where the expression for E - (|) is given by 

E - <j) = -2 tan 1 [e sin <f>/(1 + Vi - e 2 + e cos <J»)J 

4.3 TIME TERMINATION PROCEDURE 

Because the PS element set uses the true anomaly t as the independent vari- 
able, an iteration procedure is necessary to stop at a specific time t*- 
Within the ASOP program, this iteration is performed by the TIMEPS subroutine 
in the following manner: 

An expression for the derivative of time with respect to the true anomaly x 
is given in the PS theory as 


dt 

dt 


2 . 

r /q 


This expression 


can be linearly approximated by 


At 

Ax 


2 


r 


/q 


(4.3a) 
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where At = t final - t n and At = t n+1 - t n . Equation 4.3a then yields a re- 
cursive formula for refining an initial estimate of T of 

t ' = f - At q/r^ (4.3b) 

n+1 n 

Using equation 4.3b, an initial estimate of T is refined until the associated 
value of t equals the desired final time tf j_ na j_ . 

In order to start the iteration, a suitable initial value of T is necessary. 
This value is determined by first assuming that a circular orbit is being 
used. With this assumption, equation 4.2a reduces to 

% = t 


and an initial approximation for t can be written as 


T 

O 


3a 

(t-. i " 

final 4 


The value of 3a^/3T is taken from the analytical theory (ref. 5) during the 
initialization procedure . 

Therefore, the full algorithm is 

a. Set the iteration counter n to zero, and compute the initial approximation 
f “* T n = T o- 

b. Determine the PS elements at t q . 

c. Determine the time at T n . 

d. If I tfinal “ t n | < TOLerance, then STOP; otherwise n = n + 1 

e. If n > % ax , then print a diagnostic message and STOP 

f. Compute a new approximation for f n using 

T n * T n-1 " (t n-1 ‘ t final )q/r2 

g. Go to step b. 

Values for TOLerance and njjj ax have been preset within the TIMEPS subroutine 
to 10”' and 15, respectively. 
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APPENDIX A 

AVAILABLE UNITS AND PHYSICAL CONSTANTS 

Within the ASOP program, there are a number of options for the input units. 
The compatibility of the input values with the selected physical constants 
is the responsibility of the user and is controlled by the input flag IUNITS. 
The following constants are listed in order according to the value assigned 
to IUNITS. If no value for IUNITS is explicitly given, then 1 is assumed. 

R e = Earth radius (equatorial) 

= 6378.140 km (IUNITS =1,5) 

= 3443.920 ran (IUNITS =2,6) 

= 2.092566 x 10 7 ft (IUNITS = 3) 

= 6.378140 x 10 6 m (IUNITS = 4) 

= 1.0 E.r. (IUNITS = 7) 

p = Gravitational constant of the Earth 

= 3-986013 x 10 5 km 3 /sec 2 (IUNITS = 1) 

= 6.275029 x 10 1 * nm 3 /sec 2 (IUNITS = 2) 

= 1.407647 x 10 16 ft 3 /sec 2 (IUNITS = 3) 

= 3*986013 x 10 11{ m 3 /sec 2 (IUNITS = 4) 

= 5.165873 x 10 12 km 3 /hr 2 (IUNITS = 5) 

= 8.132438 x 10 11 nm 3 /hr 2 (IUNITS = 6) 

= 5.530432 x 10 -3 E.r. 3 /min 2 (IUNITS = 7) 

DAY = Time conversion 

= 8.64 x 10 4 sec/day (IUNITS = 1,2, 3, 4) 

= 1.44 x 10 3 min/day (IUNITS = 7) 

=2.40 x 10 1 hr/day (IUNITS = 5,6) 


161 



77FM50 


J 2 = J 2 coefficient of the Earth's geopotential 

= 1.082637 x 10~ 3 
e = EPS = .3/2 (.U J 2 R e 2 ) 
it = PI = 3.14159 26535 898 

If Keplerian elements are input (IEL =1), then the semimajor axis must be input 
in a distance compatible with the selected value of IUNITS, i.e., if IUNITS = 2, 
then EL(1) must be given in nautical miles; if IUNITS = 4, then EL(1) must be in 
meters, etc. All angles must be given in degrees.. 

For Cartesian coordinates, the input unit must be 


IUNITS 

EL(1) EL(3) 

EL(4) EL(6) 

1 

km 

km/ sec 

2 

nm 

nm/sec 

3 

ft 

ft /sec 

4 

m 

m/sec 

5 

km 

km/hr 

6 

nm 

nm/hr 

7 

E.r. 

E.r. /min 

All computations within 

ASOP are done using 

the input units. 
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APPENDIX B 


REQUIRED CONTROL CARDS 

1 . >§QUAL FM6-N08569 

2. >@ASG,A *NUMEG. 

3- >@XQT *NUMEG. ASOP-PROG 

4. INPUT DATA USING NAMELIST ' $INPUT' 

5 . >K$INPUT \ 

Necessary input parameters; \ 
see section 2.1.1. All l or 

parameters must be preceded 1 

by at least one space. / 

6 . >K$END 
Initial output 

7. ENTER: X = EXECUTE; S = STOP; C = 

> 

If an X is entered , the program will printout the desired information , 
and the sequence will begin again at line 4. 

If an S is entered in response to the prompt *>*, then the message 
**NORMAL PROGRAM TERMINATION** should appear. 

If a C is entered, the program will respond with **CHANGE DATA USING 
THE NAMELIST ' $INPUT’** and the sequence will begin again at line 5. 

8. >@E0F 

**NORMAL PROGRAM TERMINATION** 

> = system prompt; 
b = necessary space 


! >@ADD filename. element name 

> Input the necessary correction 

>i5$END 

CHANGE INPUT 
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APPENDIX C 

ASOP DEFAULT VALUES 

IEL = 1 (Keplerian elements) 

EL = None ; must be input by user 

STOP = 100 .ODO 

ISTOP = 2 (STOP value is in revolutions) 

PRINT = 0 .ODO. 

IPRINT = 0 (No -printout; PRINT value ignored) 

IDRAG = 1 (Drag terms desired) 

AREA = 185. 3D0 (m 2 ) j 

CD = 2.2D0 > Shuttle average 

XMASS = 90700. ODO' (kg) 1 

ILONG =1 (J 2 and short period secular terms) 

NMAX = 2 (Include J 2,0 zona ^- term) 

MMAX = 0 (No tesseral terms) 

IPSPRT = 0 (Do not print PS elements) 

IUNITS = 1 (Input and output values are given as 

distance = km 
velocity = km/ see 
time = day 
angles = deg) 
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APPENDIX D 

SUBROUTINE STORAGE REQUIREMENTS 


Subroutine 

name 

FORTRAN compiler : S0E3 

MAP processor : 28R2 RL71-3 

Storage requirements, words 

Octal 

Decimal 

MAIN 

171 

121 

AEIXYZ 

302 

194 

ASOP a 

235 

46 

CANF0R a 

524 

340 

CDTOJD 

176 

126 

COEFF 3 

322 

210 

CONST 3 

406 

262 

DENSTY 3 ’ b 

672 

442 

DETER# 

1 467 

823 

DRAG a 

115 

77 

FPRIME a 

1 177 

639 

GEOPOT 3 

430 

280 

. IL0G10 a 

106 

70 

INITAL a > b 

255 

173 

INPUT 

407 

263 

JDTOCD 

160 

112 


a ASOP subroutine package programs. 

These subroutines are models that may or may not be found in user's own 
library. 
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Subroutine 

Storage requirements, words 

name 

Octal 

Decimal 

LONGPP a 

1 705 

965 

MATiN a,b 

573 

379 

MTOECC a 

163 

115 

OUTPUT 

710 

456 

P0TEXP a 

2 214 

1 164 

PREPD a 

2 212 

1 162 

PREPS 3 }b 

276 

190 

PREPT 3 ^ 

1 031 

537 

PSANS 3 

1 712 

970 

PSTOX 3 

305 

197 

RECUR a 

163 

115 

SACT 3 » b 

545 

357 

SUN a ’ b 

153 

107 

TABLE a 

1 423 

787 

TIMEPS a 

502 

322 

TIMEXP 3 

200 

128 

XTOPS 3 

576 

382 • 

XYZAEI 

406 

262 

Subtotals 



ASOP program 

30 745 

12 773 

ASOP subroutine package 

25 747 

11 239 

COMMON block storage requirements 

13 067 

5 687 


3 ASOP subroutine package programs. 

°These subroutines are models that may or may not be found in user’s own 
library . 
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Totals 

Octal 

Decimal 

ASOP program 

44 034 

18 460 

ASOP. subroutine package 

41 036 

16 926 
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The subroutines labeled with a ’c* are models that may or may not be found in 
the user's own subroutine library. They also may be subject to change depending 
on the user's needs or requirements. For instance, PREPT initializes the coeffi- 
cients of a very accurate but extremely large geopotential model (18th order x 
18th degree). If such accuracy is not required a 4 x 4 model, for example, may 
be included to reduce storage requirements. The storage requirements for these 
standard library routines are 4211s or 2185iq words. Therefore, the storage re- 
quirements without these standard library routine options- are: 

Octal Decimal 

Totals 

ASOP program 37 623 16 275 

ASOP subroutine package 34 625 14 74 1 


Some common blocks may also be reduced if, for instance, a smaller geopotential 
model was used. This would further reduce the storage requirements. 

Values given for the storage requirements of the individual subroutines are 
values returned by the FORTRAN compiler when forming a relocatable element. The 
final, executable program will require more space because of the system library 
modules that must also be included. 

The ASOP program, with system routines and load tables, occupies over 28 000-jo 
words of storage. This has been reduced to approximately 21 500 -jq words by a 
simple overlay structure. Most of the program remains available at all times. 
Only three sections overlap each other: (1) input subroutines, (2) initi- 

alization of the drag model, and (3) computation of mean energy due to tesseral 
and sectorial geopotential harmonics. This overlay structure does not increase 
execution time significantly because the initialization routines are overlayed 
with updating routines. The core is swapped out only once unless more input i3 
brought in. a 


a For long prediction intervals (STOP 3 days) or for conditions near reentry, 
the initialization is performed more than once to account for second order 
perturbations . 
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Two test cases were run to compare the total CPU time needed to initialize 
and update an orbit prediction. Both were orbit predictions for 1 day with 
stop options on time and revolutions . The first test case consists of four 
sets of conditions: 

a. An 8th order 8th degree (8 x 8) geopotential model and a diurnal atmospheric 
drag model 

b. An 8 x 0 (no tesseral) geopotential model and a diurnal atmospheric drag 
model 

c . A 5 x 2 geopotential model and a static atmospheric drag model 

d. An 8 x 0 geopotential model with no drag 

The second test case consists of three sets of conditions . There is no drag 
because of the high eccentricity. 

a. An 8 x 8 geopotential model 

b. An 8 x 0 (no -tesseral) geopotential model 

c. A 5 x 2 geopotential model 


169 



77FM50 


Initialization is the time spent in accepting input data, transforming the 
state, and preparing to update the state. Updating is the process of actually 
propagating the state to a desired condition. Thus, if a particular problem 
requires 50 ms to initialize and 8 ms to update, then the total execution time is 
50 + 8 s 58 ms. If nine intermediate states are desired-,, then the total execu- 
tion time required to determine the nine intermediate states and one final state 
is 50 +9*8+8= 130 ms. Note that the initialization usually requires more 
time than this to update, but it is required only once a . 


a For long prediction intervals (STOP > 3 days) or for conditions near reentry, 
the initialization is performed more than once to account for second-order 
perturbations . 
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Test Case 1 : 

Small eccentricity (0.02) 

orbit 


Initial conditions: A = 6712.39 km 

W = 30° 

e = .02 
Q = 20° 

i = 30° 
M = 20° 



i 

Stop on time 

Stop on revs^ 

a. 8x8 

Initialize 

350 ms 

350 ms 


Update 

17 ms 

8 ms 

b. 8x0 

Initialize 

50 ms 

50 ms 


Update 

17 ms 

8 ms 

\ 

c. 5x2 

Initialize 

25 ms 

25 ms 


Update 

15 ms 

7 m3 

d . No drag 

Initialize 

300 ms 

300 ms 


Update 

17 ms 

8 ms 

Test Case 2: 

High eccentricity (0.72729) orbit 


Initial conditions: A = 24407.29 km 

W = 0° 

e = .72729 
ft = 20° 

i = 28.6° 
M = 0° 



Stop on time^ . 

Stop on revs 1 

a. 8x8 

Initialize 

750 ms 

750 ms 


Update 

30 ms 

30 ms 

b. 8x0 

Initialize 

18 ms 

18 ms 

• 

Update 

30 ms 

10 ms 

c. 5x2 

Initialize 

30 ms 

30 ms 


Update 

25 ms 

7 ms 

All times refer to the execution of a 
1110-EXEC 8 system. The execution time 
environment at the time of execution. 

FORTRAN V program on a UNIVAC 
will also depend on the computer 
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APPENDIX E 

GENERAL VARIABLE' ABBREVIATIONS AND DEFINITIONS 
PS elements ; 

Coordinates a 3» a 4 

(Note: a 5 = Pi i ff 6 = °7 ~ p 3> a 8 = p 4^ 

Momenta P-|, P 2 > P 3 » P 4 

Independent variable t (true anomaly) 

PS Hamiltonian: 



PS elements : 

<(> = true anomaly 
g = argument of perigee 
h = argument of the ascending node 
& = time element 

$ = conjugate to 4 >, related to the two-body energy 
G = total angular momentum 
H = Z component of the angular momentum 
L = total energy (two-body plus perturbing potential) 

(Note: For a complete description of the relationship between the DS and the PS 

elements, see reference 5.) 
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Cartesian coordinates : 

■ 

X = (X x , Xy, X z ) = position vector 
'V = (V x , V y , V z ) = velocity vector 
r = magnitude of the position vector 
t = physical time 

Keplerian elements : 

a = semimajor axis 

e = eccentricity 

i = inclination to the equator 

to = argument of peri center 

Si = argument of the ascending node 

M = mean anomaly 


Planetary variables : (see appendix A for the numerical values used) 


R e = equatorial radius 
U = gravitational constant 


General: 


km 

= 

kilometers 

min 

run 

= 

nautical miles 

rad 

ft 

= 

ft 

deg 

m 

= 

meters 

t 

E.r. 

r 

Earth radius 


sec 

n 

seconds 


hr 


hours 



minutes 

radians 

degrees 

time 

denotes a vector as X 
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APPENDIX F 

EQUATIONS OF THE ANALYTICAL THEORY 

A complete first-order solution for the motion of a satellite perturbed by 
oblateness has been developed (ref. 5). A brief outline was given in reference 17' 
and- is reproduced in this appendix. 

The Hamiltonian for the J 2 perturbed case can be written as 



and 


e = 3/2 (J 2 V R e 2 ) 


R e is the mean equatorial radius of the central body; U is the gravitational 
constant of the central body, and Jg is the Jg oblateness coefficient. 

The differential equations are solved by a method of Von-Zeipel. The elements 
undergo a canonical transformation through a determining function S-| so that 
the short periodic terms are eliminated from the Hamiltonian. The equations 
of motion in the transformed system o* may then be solved with an accuracy 
of order e . 

The solution algorithm can be divided into three steps: 
a. Canonical transformation to eliminate the short periodic terms: 





P 


k,0 


•f £ 


- e 


3S 1 

(a , 

P _ ) 

9 Pk,0 

o’ 

0 

3S, 


<v 


k,0 


k 

='1,2, 3, 4 
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b. The analytical integration of the transformed equations of motion: 

a 1 = a i,0 + V 

°2 = °P . n cos (A ? t) ' Pj n sin (A o t) 


2,0 v " 2 l/ " *^ 2,0 v ~ 2 ‘ 

a 3 = a 3,0 COS (A 3 t) ~ p 3,0 Sin (A 3 t) 

a 4 8 a 4,0 + V 
P 1 = P 1,0 

P 2 = P 2,0 ° OS ( A 2 T) + ff 2,0 Sin ( A 2 T) 

P 3 = P 3,0 cos ^ A 3 t ^ + a, 3 0 sin ^ A 3 t ^ 

P 4 8 p 4,0 


c. The back transformation: 

3S 

k k 3p k 
3S 


o, = a} - e • — 1 (CT f ,p’) 




k = 1,2, 3, 4 


If one defines 


3S 

' 1k= 55 : 


1 


where S-i = - — w y 
G2 


then 


-1 


'Ik 


/ 2 w y G \ 

(v-t y + W y k 1 ! J 


where 


w = 


Q | 

2pq 
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W k = 2 

2p q 


h [ pq Q k - Q <p k q + qp k>] 


y * JU1 <S)1 Y + T * v 

3 

y, * I i (t t \ - Yj, V * V* * r lh 

36 - 1 


= 2 <5 Ak n A + y £k C £ ) k = 2, 3,..., 8 


„ 1/2 2 . 

G = a 5 - “ (ct 2 + a 6 } 

G fc = 0 for k = 1,3, 4,7,8 

°2 = ~°2 
V 1 
°6 = ~°6 

Here p, p k » q, q k> Q, Q k> S £ , ri £ , Y £ ,. C £ and 6 £k , y £k are displayed 


1 f / 1 _ 2 _ 2\ H 

V L '\2 2 + 6 ) + ^ 


VPp 

P 2 . -2- c 2 


Vyp 

p 6 * - 2 V “6 


o Vvp 

P8 = ' (2cr 8 )3/2 
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P k = 0 for k = 1,3, 4, 5,7 


q = " o (cr 6 2 + a ? 2 ~ CT q) + — 7 = 
2 b 2 5 2^2^ 


q 2 = ~*2 


Qc = - 
5 2 


U 1 

qg " " i (2a 8 )3/2 
q = 0 for k = 1,3, 4,6, 7 



- ; (CT 2 2 - °6 2) 


x 1/2 


] 


6 2QU 2 


q = VHP 

8 2QP 2 

Ok = 0 for k = 1,3, 4, 5, 7 


6 1 = f a 6 " l (a 6 c - a 2 s) 


s a 6 *1 

6 1, = - + — B - - (0 g -as) 
12 2 3 2 2 62 22 


^ B 0 a 6 1 

5,6 ' 5 ‘ 2 T 6 - 2 ( V 6 - W 


Ik 


CT 6 1 

= r \ " 2 ( Vk " a 2 S k } for k = 1 >3, 4, 5, 7, 8 


177 



77FM50 


h ’ | °2 * \ ' V * V> 


/ B C \ cr, 1 

Y 12 = Vi *lrz* 2 * 2 ,0 « S 2 + 0 2 0 2 ) 


Y 16 = 1 * 5 S B 6 - I ( Vs * a 2°6> 


'Ik 


= -- B,, + - (a 6 s k + cx 2 c k ) for k = 1,3, 4, 5, 7, 8. 


S 2 = - -- 

2 2Q 


5 s k = 55 (| Sc-Ok)^ k = 1,2,3 8 


Yo = 


2Q 


Y 2k = " 2Q (q 5 k " s k) f "‘ k = 1,2,3 8 


{ 3 = ‘ 6 ' V * V> 


{ 32 ' ' 1 (0 2 s 2 + °6°2 + s> 


S 36 - * I <0 2 s 6 + V6 * 0) 


6 3k = " g C °2 s k + ^“k 5 for k = 


Y 3 * 6 ( V - °2 o) 
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Here 


Y 32 - g (<T 6 S 2 ' V 2 - C) 


Y 36. = g (a 6 s 6 ' °2°6 + s) 


Y 3k = g (0 6 s k " a 2°k ) for k = 1 ’ 3 ’ 1 *’ 5 ’ 7 ’ 8 


= sin £a i 


= cos £a 1 for 8, = 1,2,3 


j, s, c^, s^, B, B^, H and are displayed 


c = (G + H) 


2 2 
Orj - Og 


H,c 

c_ = 3 - (G + H) Cf 0 

3 (G + H) 3 


H 7 c 

c„ = + (G + H) cr 


7 (G + H) 


G k + H k 

c k = ~ c for k = 1 , 2 , 4 , 5 , 6, 8 

K (G + H) 


S = -(G + H) CT^ Orj 


H,s 

s, =— 3 (G + H) Orj 

3 (G + H) 7 


H 7 s 

s„ = — i (G + H) a, 

7 (G + H) - 


(G k + H k ) 

s,_ - s for k 


k (G + H) 


= 1 , 2 , 4 , 5 , 6, 8 


2 2 
B = G - 3H 
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= 2(GG^ — 3HH k ) for* k - 

H = G - - (a 2 +• a 7 2 ) 

2 3 ( 



H k = G k for k = 1,2, 4, 5, 6, 8 

Abbreviations used in the integration of the primed system 


A„ = -f„(b - -) 


y 


4 ■ 2 3 (20 8 )3/2 


A 3 = f b 3 


A 2 = 


=i [ f 2 h y^] ^ 

-/i HH 


a i ■ 1 ♦ 


f = — 
PQ 


f = -- /-yp + 2q 

1 \2 } 


f 0 = - - (W> + 2q VPp) 


f * = (2P 


L y5Tl + 2< 3 «) 
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APPENDIX G 


STANDARD 

FORTRAN VARIABLES 

USED IN ASOP 

FORTRAN 

variable 

Program 

location 

Description 

AEI ( 6 ) 

OUTPUT 

Character array to accompany Keplerian element 
output 

ANG(3) 

OUTPUT 

Character array to accompany any angular output 

B(3,3) 

AEIXYZ 

Keplerian elements to Cartesian coordinates 
transformation matrix 

BS 

/S 1 STAV/ 

1 - H/G 

BY3 

/CONSTW/ 

1/3 

BY 6 

/CONSTW/ 

1/6 

C( 8 ) 

PSANS 

9c/3ok > Sc/Spjc , k = 1 ,2,3,4 

CHECK 

OUTPUT 

Energy check value 

CINC 

AEIXYZ 

Cosine of the orbital inclination with respect 
to the Earth's equator (cos i) 

CN 

/CONSTW/ 

+1 depending on value of NN 

CNODE 

AEIXYZ 

Cosine of the argument of the ascending node 
(cos ft) 

COMEGA 

AEIXYZ 

Cosine of the argument of pericenter (cos w) 

COSEA 

AEIXYZ 

Cosine of the eccentric anomaly (cos E) 

C0SFC2 

PSANS 

cos (A 2 t) (see appendix F) 

COSFC3 

PSANS 

cos (A 3 T) (see appendix F) 

CS 

PSANS 

( 1 / 2 ) (G + H) (P 3 2 + 0 3 2) = small »c' 

DAYS(7) 

CONST 

Storage array of possible values of DAY 

DAYS 

OUTPUT 

Print value (DAYS = TIME/DAY) 

DAYS 

TIMEPS 

Total days elapsed 
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FORTRAN 

variable 

Program 

location 

Description 

DELTA 1 

PSANS 

«i 

DELTA2 

PSANS 

«2 

DELTA3 

PSANS 

63 

DELI ( 8 ) 

PSANS 

36.,/aak, ss^apk k = 1 , 2 , 3 , 4 

DEL2(8) 

PSANS 

862 / 90 ^’ 982/3Pk k = 1 ,2,3,4 

DEL3(8) 

PSANS 

363 / 80 ^, 383/3pk k = 1,2,3 ,4 

DSB(4) 

PSANS 

3b/38 k k = 1,2,3, 4 

DSF(4) 

PSANS 

3f/38 k k = 1 , 2 , 3, 4 

DST ( 7 ) 

OUTPUT 

Character array to accompany any distance 
output 

EA 

/AEIXYZ) 

(XYZAEIJ 

The eccentric anomaly of the satellite computed 
from Kepler's equation (E) (rad) 

EAO 

AEIXYZ 

Old value of EA; used when iterating to solve 
Kepler's equation (rad) 

ECOSE 

XYZAEI 

E cos e 

ECOSPH 

( PSTOX \ 
\XTOPS j 

E cos <j> 

EL(6) 

INPUT 

Initial conditions of the satellite given in 
Keplerian elements or Cartesian coordinates; on 
output , it will contain the Keplerian -elements . 
EL( 1 ) X or a 

(2) Y or e 

(3) Z or i 

(4) X or to 

(5) I or S! 

(6) Z or M 

EMINPH 

f PSTOX \ 
\XTOPS j 

E - $ 

EROOT 

/ PSANS3/ 

V 2 P4 P/P 

EROOT 

XTOPS 

Vi . - 2Q($ - G) 
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FORTRAN 

variable 

Program 

location 

Description 

ESINE 

XYZAEI 

E sin e 

ESINPH 

fPSTOX\ 

\XTOPS) 

E sin <j> 

ETA2 

PSANS 

sin 2a i 

ETA3 

PSANS . 

sin 3<?-| 

FACTOR (4) 

PSANS 

Derivatives of the DS Hamiltonian and its 
combinations (A-| , A 2 , Ag, Aj|) (see appendix F) 

FN(19,19) 

POTEXP 

Inclination function 

FS 

/S1STAV/ 

f 

FSSQ 

PSANS 

f 2 

GAM1A1 

PSANS 

Y 1 

GAMMA2 

PSANS 

Y 2 

GAMMA3 

PSANS 


GAM1(8) 

PSANS 

3y-|/3aj^, Sy^/Spjj k = 1 ,2,3,4 

GAM2( 8 ) 

PSANS 

3y 2 /3o’ k , 3y 2 /3p k k = 1,2, 3, 4 

GAM3(8) 

PSANS 

3y 3 /3a k , 3y 3 /3p k k = 1,2, 3,4 

GC(8) 

/S1STAD/ 

3G/3a k ,‘ 3G/3p k k = 1,2, 3, 4 

GCAP( 17 ) 

POTEXP 

Eccentricity function 

GCIN 

PSTOX 

G” 1 

GCSQ 

XTOPS 

G 2 

GIN 

/S1STAV/ 

G -1 

G1NSQ 

/S1STAV/ 

G" 2 

GM3H 

PSANS 

G - 3H * 

GPH ' 

/S1STAV/ 

G + H 

GSQ 

PSANS 

G 2 
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FORTRAN 

variable 

Program 

location 

Description 

G1 

XTOPS 

YV Z - ZV y = G x 

G1SQ 

XTOPS 

r 2 
G x 

G2 

XTOPS 

ZV X - XV Z = G y 

G2SQ 

XTOPS 


G3 

XTOPS 

XVy - = G z 

H 

XXZAEI 

Total angular momentum 

HC(8) 

/S1STAD/ 

3H/3a k , 3H/3p k k = 1,2, 3, 4 

HMS(li) 

OUTPUT 

Character array to accompany any time output 

HOG 

/SISTAV/ 

H/G 

HSQ 

PSANS 

H 2 

IA(3) 

OUTPUT 

Character array of blanks and asterisks 

I EL 

INPUT 

Flag to determine if input values of EL are given 
as Keplerian elements or Cartesian coordinates 

= 1 Keplerian 

= 2 Cartesian 

XERR 

TIMEPS 

Error counter 

IFORM 

OUTPUT 

Flag to determine if initial or final condition 
messages is to be printed 

= 1 initial condition message 

= 2 no message (intermediate print) 

= 3 final condition message 

IMARK 

PSANS 

Flag determining if one or two passes have been 
made 

=1 1st pass 
= 2 2nd pass 
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fortran 

variable 

Program 

location Description 

ITER 

TIMEPS 

Total number of iterations allowed 

IXP 

OUTPUT 

Pointer to the IA array 



IFORM i 3 ♦ IXP = IP 



IFORM = 3 + IXP = ISTOP 

L 

PSANS | 


LC 

XTOPS ) 

ii 

<? 

it 

LS 

PSANS 

& = OH 

NEWX j 

f MAIN \ 
[ ASOP / 

Flag to determine if ASOP program is to be 
initialized 



= 0 no 



= 1 yes 

NN 

PSANS 

Flag to determine if initializing or computing 



= 0 initializing 



= 1 computing 

P 

XYZAEI 

1 - e 2 

P(8) 

/S1STAD/ 3p/9c^, 3p/3pic k = 1,2, 3, 4. 

PHI 

psans ; 

| 

PHIC 

XTOPS . 

> $ = Pi = o 5 

1 

f GEOPOT' 

\ 

POT •< 

i 

PSTOX 
1 XTOPS 

/ Magnitude of Earth's gravitational potential 
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FORTRAN 

variable 

Program 

location 

Description 

p 5 I 

[PSANS \ 
[PSTOX ) 

1/2 

Q(8) 

/S1STAD/ 

3q/ac k , 3q/ap k k = 1,2, 3, 4 

QCIN 

PSANS 

Q- 1 ' 

QCSQ 

PSANS 

Q 2 

QCV(8) 

/S1STAD/ 

3Q/3c k , 3Q/ap k k = 1,2, 3, 4 

QS 

/PSANS2/ 

1/2 (y- . g ) 

\2L / 

RCAP 

f PSTOX \ 

R (see section 4.2) 


Utops j 


RCAPDT 

PSTOX 

XTOPS 

3R/3t (see section 4.2) 

RCOSF 

XYZAEI 

r cos f; f = true anomaly, r = R 

RCOSL 

XYZAEI 

r cos L; L = mean anomaly, r = R 

RDOT 

PSTOX 

Magnitude of velocity vector 

RES( 7 ) 

CONST 

Storage array of possible values of RE 

REVS 

OUTPUT 

Total number of revolutions predicted 
(REVS = TAU/2 tt) 

ROP 

PSTOX 

r/p (see section 4.2) 

RRDOT 

XYZAEI 

X * V ‘ 

RSINF 

XYZAEI 

\ 

r sin f, f = true anomaly, r = R 

RSINL 

XYZAEI 

r sin L, L = mean anomaly, r = R 

RSQ 

XTOPS 

R 2 

R2I 

GEOPOT 

1/R 2 

S(8) 

PSANS • 

as/3cj k , 3s/ap k k = 1,2, 3, 4 

SINC 

AEIXYZ 

Sine of the orbital inclination with respect 
to the Earth's equator (sin i) 
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fortran 

variable 

Program 

location 

Description 

SIGINI(8) 

PSANS 

Initial values of o's and p's 

SINEA 

AEIXYZ 

Sine of the eccentric anomaly (sin E) 

SINFC2 

PSANS 

sin (A2 t) (see appendix F) 

SINFC3 

PSANS 

sin (Ag t) (see appendix F) 

SINWF 

XYZAEI 

sin (w + f) to = argument of perigee, 
f = true anomaly 

SNODE 

AEIXYZ 

Sine of the argument of the ascending node 
( sin &) 

SOMEGA 

AEIXYZ 

Sine of the argument of peri center (sin to) 

SQTGHI 

XTOPS 

-V2/(G + H) 

SS 

PSANS 

’ s 

STOPDT 

MAIN 

Value at .which next intermediate printout is 
desired (needed only if IPRINT >1) 

SUM2 

/PSANS2/ 

1/2 (<? 2 2 + P 2 2 ) = $ - G 

SUM3 

/PSANS2/ 

1/2 (a 3 2 + P3 2 ) = 2 (G - H) 

SI (4) 

PSANS 

Derivatives of the generating function' S 1 

TFIN 

TIMEPS 

Final time desired for stopping the iteration 

TW03 

/CONSTW/ 

2/3 

VEL ( 7 ) 

OUTPUT 

Character array to accompany any velocity 
output 

VI (3) 

AEIXYZ 

Velocity vector with respect to the orbital 
plane (V x » , V y ' , V z * = 0) 

VSQ 

XYZAEI 

Magnitude of the velocity vector, squared 
(V 2 ) 

W(8) 

PSANS 

9w/3a k , 9w/3p k k = 1,2, 3 , 4 

WS 

PSANS 

w 
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FORTRAN 

Program 


variable 

location 

Description 

WX j 


Components of the total angular momentum 

WY > 

XYZAEI 

(WX, WY, WZ) (distance 2 /time) 

WZ J 



XMUS(7) 

CONST 

Storage array of possible values of XMU 

XIN(8) 

ASOP 

Identical to X but allows ASOP subroutine to 
be removed from the stand alone program 

XYZ(6) 

OUTPUT . 

Character array to accompany output of the 
Cartesian state vector 

X1(3) 

AEIXYZ 

Position vector of the satellite with respect 
to the orbital plane (X 1 , Y', Z' = 0) 

X3ROOT 

/PSANS3/ 

\yUG - Pg 2 - pg 2 j /G (see section 4.2) 

Y 

PSANS 

3y/3a k ., 3y/3p k k = 1,2, 3, 4 

YS 

PSANS 

y 

ZCAP1 

XTOPS 

Z 1 

ZCAP2 

XTOPS 

z 2 

ZET2 

PSANS 

cos 2<J-j 

ZET3 

PSANS 

cos 3a.| 
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Figure 41.- Data flow in ASOP-general subroutine linkage. 


a GEOPOT subroutines are the same 
^ILOGIO subroutines are the same 
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Figure 42.- General subroutine linkage in 
removable ASOP subroutine package. 


a GE0P0T subroutines are the same 
^ILOGIO subroutines are the same 
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